chiark / gitweb /
version, copyright-year bump
[nlopt.git] / api / optimize.c
1 /* Copyright (c) 2007-2010 Massachusetts Institute of Technology
2  *
3  * Permission is hereby granted, free of charge, to any person obtaining
4  * a copy of this software and associated documentation files (the
5  * "Software"), to deal in the Software without restriction, including
6  * without limitation the rights to use, copy, modify, merge, publish,
7  * distribute, sublicense, and/or sell copies of the Software, and to
8  * permit persons to whom the Software is furnished to do so, subject to
9  * the following conditions:
10  * 
11  * The above copyright notice and this permission notice shall be
12  * included in all copies or substantial portions of the Software.
13  * 
14  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
15  * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
16  * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
17  * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
18  * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
19  * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
20  * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. 
21  */
22
23 #include <stdlib.h>
24 #include <math.h>
25 #include <float.h>
26
27 #include "nlopt-internal.h"
28
29 /*********************************************************************/
30
31 #ifndef HAVE_ISNAN
32 static int my_isnan(double x) { return x != x; }
33 #  define isnan my_isnan
34 #endif
35
36 /*********************************************************************/
37
38 #include "praxis.h"
39 #include "direct.h"
40
41 #ifdef WITH_CXX
42 #  include "stogo.h"
43 #endif
44
45 #include "cdirect.h"
46
47 #ifdef WITH_NOCEDAL
48 #  include "l-bfgs-b.h"
49 #endif
50
51 #include "luksan.h"
52
53 #include "crs.h"
54
55 #include "mlsl.h"
56 #include "mma.h"
57 #include "cobyla.h"
58 #include "newuoa.h"
59 #include "neldermead.h"
60 #include "auglag.h"
61 #include "bobyqa.h"
62 #include "isres.h"
63
64 /*********************************************************************/
65
66 typedef struct {
67      nlopt_func f;
68      void *f_data;
69      const double *lb, *ub;
70 } nlopt_data;
71
72 #include "praxis.h"
73
74 static double f_bound(int n, const double *x, void *data_)
75 {
76      int i;
77      nlopt_data *data = (nlopt_data *) data_;
78      double f;
79
80      /* some methods do not support bound constraints, but support
81         discontinuous objectives so we can just return Inf for invalid x */
82      for (i = 0; i < n; ++i)
83           if (x[i] < data->lb[i] || x[i] > data->ub[i])
84                return HUGE_VAL;
85
86      f = data->f((unsigned) n, x, NULL, data->f_data);
87      return (isnan(f) || nlopt_isinf(f) ? HUGE_VAL : f);
88 }
89
90 static double f_noderiv(int n, const double *x, void *data_)
91 {
92      nlopt_data *data = (nlopt_data *) data_;
93      return data->f((unsigned) n, x, NULL, data->f_data);
94 }
95
96 static double f_direct(int n, const double *x, int *undefined, void *data_)
97 {
98      nlopt_data *data = (nlopt_data *) data_;
99      double f;
100      f = data->f((unsigned) n, x, NULL, data->f_data);
101      *undefined = isnan(f) || nlopt_isinf(f);
102      return f;
103 }
104
105 /*********************************************************************/
106
107 /* get min(dx) for algorithms requiring a scalar initial step size */
108 static nlopt_result initial_step(nlopt_opt opt, const double *x, double *step)
109 {
110      unsigned freedx = 0, i;
111
112      if (!opt->dx) {
113           freedx = 1;
114           if (nlopt_set_default_initial_step(opt, x) != NLOPT_SUCCESS)
115                return NLOPT_OUT_OF_MEMORY;
116      }
117
118      *step = HUGE_VAL;
119      for (i = 0; i < opt->n; ++i)
120           if (*step > fabs(opt->dx[i]))
121                *step = fabs(opt->dx[i]);
122
123      if (freedx) { free(opt->dx); opt->dx = NULL; }
124      return NLOPT_SUCCESS;
125 }
126
127 /*********************************************************************/
128
129 /* return true if [lb,ub] is finite in every dimension (n dimensions) */
130 static int finite_domain(unsigned n, const double *lb, const double *ub)
131 {
132      unsigned i;
133      for (i = 0; i < n; ++i)
134           if (nlopt_isinf(ub[i] - lb[i])) return 0;
135      return 1;
136 }
137
138 /*********************************************************************/
139
140 #define POP(defaultpop) (opt->stochastic_population > 0 ?               \
141                          opt->stochastic_population :                   \
142                          (nlopt_stochastic_population > 0 ?             \
143                           nlopt_stochastic_population : (defaultpop)))
144
145 nlopt_result nlopt_optimize(nlopt_opt opt, double *x, double *minf)
146 {
147      const double *lb, *ub;
148      nlopt_algorithm algorithm;
149      nlopt_func f; void *f_data;
150      unsigned n, i;
151      int ni;
152      nlopt_data d;
153      nlopt_stopping stop;
154
155      if (!opt || !x || !minf || !opt->f) return NLOPT_INVALID_ARGS;
156
157      /* copy a few params to local vars for convenience */
158      n = opt->n;
159      ni = (int) n; /* most of the subroutines take "int" arg */
160      lb = opt->lb; ub = opt->ub;
161      algorithm = opt->algorithm;
162      f = opt->f; f_data = opt->f_data;
163
164      if (n == 0) { /* trivial case: no degrees of freedom */
165           *minf = f(n, x, NULL, f_data);
166           return NLOPT_SUCCESS;
167      }
168
169      *minf = HUGE_VAL;
170      
171      d.f = f;
172      d.f_data = f_data;
173      d.lb = lb;
174      d.ub = ub;
175
176      /* make sure rand generator is inited */
177      if (!nlopt_srand_called)
178           nlopt_srand_time(); /* default is non-deterministic */
179
180      /* check bound constraints */
181      for (i = 0; i < n; ++i)
182           if (lb[i] > ub[i] || x[i] < lb[i] || x[i] > ub[i])
183                return NLOPT_INVALID_ARGS;
184
185      stop.n = n;
186      stop.minf_max = opt->minf_max;
187      stop.ftol_rel = opt->ftol_rel;
188      stop.ftol_abs = opt->ftol_abs;
189      stop.xtol_rel = opt->xtol_rel;
190      stop.xtol_abs = opt->xtol_abs;
191      stop.nevals = 0;
192      stop.maxeval = opt->maxeval;
193      stop.maxtime = opt->maxtime;
194      stop.start = nlopt_seconds();
195
196      switch (algorithm) {
197          case NLOPT_GN_DIRECT:
198          case NLOPT_GN_DIRECT_L: 
199          case NLOPT_GN_DIRECT_L_RAND: 
200               if (!finite_domain(n, lb, ub)) return NLOPT_INVALID_ARGS;
201               return cdirect(ni, f, f_data, 
202                              lb, ub, x, minf, &stop, 0.0, 
203                              (algorithm != NLOPT_GN_DIRECT)
204                              + 3 * (algorithm == NLOPT_GN_DIRECT_L_RAND 
205                                     ? 2 : (algorithm != NLOPT_GN_DIRECT))
206                              + 9 * (algorithm == NLOPT_GN_DIRECT_L_RAND 
207                                     ? 1 : (algorithm != NLOPT_GN_DIRECT)));
208               
209          case NLOPT_GN_DIRECT_NOSCAL:
210          case NLOPT_GN_DIRECT_L_NOSCAL: 
211          case NLOPT_GN_DIRECT_L_RAND_NOSCAL: 
212               if (!finite_domain(n, lb, ub)) return NLOPT_INVALID_ARGS;
213               return cdirect_unscaled(ni, f, f_data, lb, ub, x, minf, 
214                                       &stop, 0.0, 
215                                       (algorithm != NLOPT_GN_DIRECT)
216                                       + 3 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 2 : (algorithm != NLOPT_GN_DIRECT))
217                                       + 9 * (algorithm == NLOPT_GN_DIRECT_L_RAND ? 1 : (algorithm != NLOPT_GN_DIRECT)));
218               
219          case NLOPT_GN_ORIG_DIRECT:
220          case NLOPT_GN_ORIG_DIRECT_L: 
221               if (!finite_domain(n, lb, ub)) return NLOPT_INVALID_ARGS;
222               switch (direct_optimize(f_direct, &d, ni, lb, ub, x, minf,
223                                       stop.maxeval, -1, 0.0, 0.0,
224                                       pow(stop.xtol_rel, (double) n), -1.0,
225                                       stop.minf_max, 0.0,
226                                       NULL, 
227                                       algorithm == NLOPT_GN_ORIG_DIRECT
228                                       ? DIRECT_ORIGINAL
229                                       : DIRECT_GABLONSKY)) {
230                   case DIRECT_INVALID_BOUNDS:
231                   case DIRECT_MAXFEVAL_TOOBIG:
232                   case DIRECT_INVALID_ARGS:
233                        return NLOPT_INVALID_ARGS;
234                   case DIRECT_INIT_FAILED:
235                   case DIRECT_SAMPLEPOINTS_FAILED:
236                   case DIRECT_SAMPLE_FAILED:
237                        return NLOPT_FAILURE;
238                   case DIRECT_MAXFEVAL_EXCEEDED:
239                   case DIRECT_MAXITER_EXCEEDED:
240                        return NLOPT_MAXEVAL_REACHED;
241                   case DIRECT_GLOBAL_FOUND:
242                        return NLOPT_MINF_MAX_REACHED;
243                   case DIRECT_VOLTOL:
244                   case DIRECT_SIGMATOL:
245                        return NLOPT_XTOL_REACHED;
246                   case DIRECT_OUT_OF_MEMORY:
247                        return NLOPT_OUT_OF_MEMORY;
248               break;
249          }
250
251          case NLOPT_GD_STOGO:
252          case NLOPT_GD_STOGO_RAND:
253 #ifdef WITH_CXX
254               if (!finite_domain(n, lb, ub)) return NLOPT_INVALID_ARGS;
255               if (!stogo_minimize(ni, f, f_data, x, minf, lb, ub, &stop,
256                                   algorithm == NLOPT_GD_STOGO
257                                   ? 0 : (int) POP(2*n)))
258                    return NLOPT_FAILURE;
259               break;
260 #else
261               return NLOPT_FAILURE;
262 #endif
263
264 #if 0
265               /* lacking a free/open-source license, we no longer use
266                  Rowan's code, and instead use by "sbplx" re-implementation */
267          case NLOPT_LN_SUBPLEX: {
268               int iret, freedx = 0;
269               if (!opt->dx) {
270                    freedx = 1;
271                    if (nlopt_set_default_initial_step(opt, x) != NLOPT_SUCCESS)
272                         return NLOPT_OUT_OF_MEMORY;
273               }                
274               iret = nlopt_subplex(f_bound, minf, x, n, &d, &stop, opt->dx);
275               if (freedx) { free(opt->dx); opt->dx = NULL; }
276               switch (iret) {
277                   case -2: return NLOPT_INVALID_ARGS;
278                   case -10: return NLOPT_MAXTIME_REACHED;
279                   case -1: return NLOPT_MAXEVAL_REACHED;
280                   case 0: return NLOPT_XTOL_REACHED;
281                   case 1: return NLOPT_SUCCESS;
282                   case 2: return NLOPT_MINF_MAX_REACHED;
283                   case 20: return NLOPT_FTOL_REACHED;
284                   case -200: return NLOPT_OUT_OF_MEMORY;
285                   default: return NLOPT_FAILURE; /* unknown return code */
286               }
287               break;
288          }
289 #endif
290
291          case NLOPT_LN_PRAXIS: {
292               double step;
293               if (initial_step(opt, x, &step) != NLOPT_SUCCESS)
294                    return NLOPT_OUT_OF_MEMORY;
295               return praxis_(0.0, DBL_EPSILON, 
296                              step, ni, x, f_bound, &d, &stop, minf);
297          }
298
299 #ifdef WITH_NOCEDAL
300          case NLOPT_LD_LBFGS_NOCEDAL: {
301               int iret, *nbd = (int *) malloc(sizeof(int) * n);
302               if (!nbd) return NLOPT_OUT_OF_MEMORY;
303               for (i = 0; i < n; ++i) {
304                    int linf = nlopt_isinf(lb[i]) && lb[i] < 0;
305                    int uinf = nlopt_isinf(ub[i]) && ub[i] > 0;
306                    nbd[i] = linf && uinf ? 0 : (uinf ? 1 : (linf ? 3 : 2));
307               }
308               iret = lbfgsb_minimize(ni, f, f_data, x, nbd, lb, ub,
309                                      MIN(ni, 5), 0.0, stop.ftol_rel, 
310                                      stop.xtol_abs[0] > 0 ? stop.xtol_abs[0]
311                                      : stop.xtol_rel,
312                                      stop.maxeval);
313               free(nbd);
314               if (iret <= 0) {
315                    switch (iret) {
316                        case -1: return NLOPT_INVALID_ARGS;
317                        case -2: default: return NLOPT_FAILURE;
318                    }
319               }
320               else {
321                    *minf = f(n, x, NULL, f_data);
322                    switch (iret) {
323                        case 5: return NLOPT_MAXEVAL_REACHED;
324                        case 2: return NLOPT_XTOL_REACHED;
325                        case 1: return NLOPT_FTOL_REACHED;
326                        default: return NLOPT_SUCCESS;
327                    }
328               }
329               break;
330          }
331 #endif
332
333          case NLOPT_LD_LBFGS: 
334               return luksan_plis(ni, f, f_data, lb, ub, x, minf, &stop);
335
336          case NLOPT_LD_VAR1: 
337          case NLOPT_LD_VAR2: 
338               return luksan_plip(ni, f, f_data, lb, ub, x, minf, &stop,
339                    algorithm == NLOPT_LD_VAR1 ? 1 : 2);
340
341          case NLOPT_LD_TNEWTON: 
342          case NLOPT_LD_TNEWTON_RESTART: 
343          case NLOPT_LD_TNEWTON_PRECOND: 
344          case NLOPT_LD_TNEWTON_PRECOND_RESTART: 
345               return luksan_pnet(ni, f, f_data, lb, ub, x, minf, &stop,
346                                  1 + (algorithm - NLOPT_LD_TNEWTON) % 2,
347                                  1 + (algorithm - NLOPT_LD_TNEWTON) / 2);
348
349          case NLOPT_GN_CRS2_LM:
350               if (!finite_domain(n, lb, ub)) return NLOPT_INVALID_ARGS;
351               return crs_minimize(ni, f, f_data, lb, ub, x, minf, &stop, 
352                                   (int) POP(0), 0);
353
354          case NLOPT_GN_MLSL:
355          case NLOPT_GD_MLSL:
356          case NLOPT_GN_MLSL_LDS:
357          case NLOPT_GD_MLSL_LDS: {
358               if (!finite_domain(n, lb, ub)) return NLOPT_INVALID_ARGS;
359               nlopt_opt local_opt = opt->local_opt;
360               nlopt_result ret;
361               if (!local_opt) { /* default */
362                    nlopt_algorithm local_alg = (algorithm == NLOPT_GN_MLSL ||
363                                                 algorithm == NLOPT_GN_MLSL_LDS)
364                         ? nlopt_local_search_alg_nonderiv
365                         : nlopt_local_search_alg_deriv;
366                    /* don't call MLSL recursively! */
367                    if (local_alg >= NLOPT_GN_MLSL
368                        && local_alg <= NLOPT_GD_MLSL_LDS)
369                         local_alg = (algorithm == NLOPT_GN_MLSL ||
370                                      algorithm == NLOPT_GN_MLSL_LDS)
371                              ? NLOPT_LN_COBYLA : NLOPT_LD_MMA;
372                    local_opt = nlopt_create(local_alg, n);
373                    if (!local_opt) return NLOPT_FAILURE;
374                    nlopt_set_ftol_rel(local_opt, opt->ftol_rel);
375                    nlopt_set_ftol_abs(local_opt, opt->ftol_abs);
376                    nlopt_set_xtol_rel(local_opt, opt->xtol_rel);
377                    nlopt_set_xtol_abs(local_opt, opt->xtol_abs);
378                    nlopt_set_maxeval(local_opt, nlopt_local_search_maxeval);
379                    nlopt_set_initial_step(local_opt, opt->dx);
380               }
381               for (i = 0; i < n && stop.xtol_abs[i] > 0; ++i) ;
382               if (local_opt->ftol_rel <= 0 && local_opt->ftol_abs <= 0 &&
383                   local_opt->xtol_rel <= 0 && i < n) {
384                    /* it is not sensible to call MLSL without *some*
385                       nonzero tolerance for the local search */
386                    nlopt_set_ftol_rel(local_opt, 1e-15);
387                    nlopt_set_xtol_rel(local_opt, 1e-7);
388               }
389               ret = mlsl_minimize(ni, f, f_data, lb, ub, x, minf, &stop,
390                                   local_opt, (int) POP(0),
391                                   algorithm >= NLOPT_GN_MLSL_LDS);
392               if (!opt->local_opt) nlopt_destroy(local_opt);
393               return ret;
394          }
395
396          case NLOPT_LD_MMA: {
397               nlopt_opt dual_opt;
398               nlopt_result ret;
399 #define LO(param, def) (opt->local_opt ? opt->local_opt->param : (def))
400               dual_opt = nlopt_create(LO(algorithm,
401                                          nlopt_local_search_alg_deriv),
402                                       opt->m);
403               if (!dual_opt) return NLOPT_FAILURE;
404               nlopt_set_ftol_rel(dual_opt, LO(ftol_rel, 1e-12));
405               nlopt_set_ftol_abs(dual_opt, LO(ftol_abs, 0.0));
406               nlopt_set_maxeval(dual_opt, LO(maxeval, 100000));
407 #undef LO
408
409               ret = mma_minimize(ni, f, f_data, (int) (opt->m), opt->fc,
410                                  lb, ub, x, minf, &stop, dual_opt);
411               nlopt_destroy(dual_opt);
412               return ret;
413          }
414
415          case NLOPT_LN_COBYLA: {
416               double step;
417               if (initial_step(opt, x, &step) != NLOPT_SUCCESS)
418                    return NLOPT_OUT_OF_MEMORY;
419               return cobyla_minimize(ni, f, f_data, 
420                                      opt->m, opt->fc,
421                                      lb, ub, x, minf, &stop,
422                                      step);
423          }
424                                      
425          case NLOPT_LN_NEWUOA: {
426               double step;
427               if (initial_step(opt, x, &step) != NLOPT_SUCCESS)
428                    return NLOPT_OUT_OF_MEMORY;
429               return newuoa(ni, 2*n+1, x, 0, 0, step,
430                             &stop, minf, f_noderiv, &d);
431          }
432                                      
433          case NLOPT_LN_NEWUOA_BOUND: {
434               double step;
435               if (initial_step(opt, x, &step) != NLOPT_SUCCESS)
436                    return NLOPT_OUT_OF_MEMORY;
437               return newuoa(ni, 2*n+1, x, lb, ub, step,
438                             &stop, minf, f_noderiv, &d);
439          }
440
441          case NLOPT_LN_BOBYQA: {
442               double step;
443               if (initial_step(opt, x, &step) != NLOPT_SUCCESS)
444                    return NLOPT_OUT_OF_MEMORY;
445               return bobyqa(ni, 2*n+1, x, lb, ub, step,
446                             &stop, minf, f_noderiv, &d);
447          }
448
449          case NLOPT_LN_NELDERMEAD: 
450          case NLOPT_LN_SBPLX: 
451          {
452               nlopt_result ret;
453               int freedx = 0;
454               if (!opt->dx) {
455                    freedx = 1;
456                    if (nlopt_set_default_initial_step(opt, x) != NLOPT_SUCCESS)
457                         return NLOPT_OUT_OF_MEMORY;
458               }
459               if (algorithm == NLOPT_LN_NELDERMEAD)
460                    ret= nldrmd_minimize(ni,f,f_data,lb,ub,x,minf,opt->dx,&stop);
461               else
462                    ret= sbplx_minimize(ni,f,f_data,lb,ub,x,minf,opt->dx,&stop);
463               if (freedx) { free(opt->dx); opt->dx = NULL; }
464               return ret;
465          }
466
467          case NLOPT_LN_AUGLAG:
468          case NLOPT_LN_AUGLAG_EQ:
469          case NLOPT_LD_AUGLAG:
470          case NLOPT_LD_AUGLAG_EQ: {
471               nlopt_opt local_opt = opt->local_opt;
472               nlopt_result ret;
473               if (!local_opt) { /* default */
474                    local_opt = nlopt_create(
475                         algorithm == NLOPT_LN_AUGLAG || 
476                         algorithm == NLOPT_LN_AUGLAG_EQ
477                         ? nlopt_local_search_alg_nonderiv
478                         : nlopt_local_search_alg_deriv, n);
479                    if (!local_opt) return NLOPT_FAILURE;
480                    nlopt_set_ftol_rel(local_opt, opt->ftol_rel);
481                    nlopt_set_ftol_abs(local_opt, opt->ftol_abs);
482                    nlopt_set_xtol_rel(local_opt, opt->xtol_rel);
483                    nlopt_set_xtol_abs(local_opt, opt->xtol_abs);
484                    nlopt_set_maxeval(local_opt, nlopt_local_search_maxeval);
485                    nlopt_set_initial_step(local_opt, opt->dx);
486               }
487               ret = auglag_minimize(ni, f, f_data, 
488                                     opt->m, opt->fc, 
489                                     opt->p, opt->h,
490                                     lb, ub, x, minf, &stop,
491                                     local_opt,
492                                     algorithm == NLOPT_LN_AUGLAG_EQ
493                                     || algorithm == NLOPT_LD_AUGLAG_EQ);
494               if (!opt->local_opt) nlopt_destroy(local_opt);
495               return ret;
496          }
497
498          case NLOPT_GN_ISRES:
499               if (!finite_domain(n, lb, ub)) return NLOPT_INVALID_ARGS;
500               return isres_minimize(ni, f, f_data, 
501                                     (int) (opt->m), opt->fc,
502                                     (int) (opt->p), opt->h,
503                                     lb, ub, x, minf, &stop,
504                                     (int) POP(0));
505
506          default:
507               return NLOPT_INVALID_ARGS;
508      }
509
510      return NLOPT_SUCCESS; /* never reached */
511 }
512
513 /*********************************************************************/
514
515 nlopt_result nlopt_optimize_limited(nlopt_opt opt, double *x, double *minf,
516                                     int maxeval, double maxtime)
517 {
518      int save_maxeval;
519      double save_maxtime;
520      nlopt_result ret;
521
522      if (!opt) return NLOPT_INVALID_ARGS;
523
524      save_maxeval = nlopt_get_maxeval(opt);
525      save_maxtime = nlopt_get_maxtime(opt);
526      
527      /* override opt limits if maxeval and/or maxtime are more stringent */
528      if (save_maxeval <= 0 || (maxeval > 0 && maxeval < save_maxeval))
529           nlopt_set_maxeval(opt, maxeval);
530      if (save_maxtime <= 0 || (maxtime > 0 && maxtime < save_maxtime))
531           nlopt_set_maxtime(opt, maxtime);
532
533      ret = nlopt_optimize(opt, x, minf);
534
535      nlopt_set_maxeval(opt, save_maxeval);
536      nlopt_set_maxtime(opt, save_maxtime);
537
538      return ret;
539 }
540
541 /*********************************************************************/