chiark / gitweb /
added CCSAQ algorithm; internal support for preconditioners (untested) not yet export...
[nlopt.git] / api / options.c
1 /* Copyright (c) 2007-2011 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 <stdio.h>
24 #include <stdlib.h>
25 #include <math.h>
26 #include <string.h>
27 #include <float.h>
28
29 #include "nlopt-internal.h"
30
31 /*************************************************************************/
32
33 void NLOPT_STDCALL nlopt_destroy(nlopt_opt opt)
34 {
35      if (opt) {
36           unsigned i;
37           if (opt->munge_on_destroy) {
38                nlopt_munge munge = opt->munge_on_destroy;
39                munge(opt->f_data);
40                for (i = 0; i < opt->m; ++i)
41                     munge(opt->fc[i].f_data);
42                for (i = 0; i < opt->p; ++i)
43                     munge(opt->h[i].f_data);
44           }
45           for (i = 0; i < opt->m; ++i)
46                free(opt->fc[i].tol);
47           for (i = 0; i < opt->p; ++i)
48                free(opt->h[i].tol);
49           free(opt->lb); free(opt->ub);
50           free(opt->xtol_abs);
51           free(opt->fc);
52           free(opt->h);
53           nlopt_destroy(opt->local_opt);
54           free(opt->dx);
55           free(opt->work);
56           free(opt);
57      }
58 }
59
60 nlopt_opt NLOPT_STDCALL nlopt_create(nlopt_algorithm algorithm, unsigned n)
61 {
62      nlopt_opt opt;
63
64      if (((int) algorithm) < 0 || algorithm >= NLOPT_NUM_ALGORITHMS)
65           return NULL;
66
67      opt = (nlopt_opt) malloc(sizeof(struct nlopt_opt_s));
68      if (opt) {
69           opt->algorithm = algorithm;
70           opt->n = n;
71           opt->f = NULL; opt->f_data = NULL;
72           opt->maximize = 0;
73           opt->munge_on_destroy = opt->munge_on_copy = NULL;
74
75           opt->lb = opt->ub = NULL;
76           opt->m = opt->m_alloc = 0;
77           opt->fc = NULL;
78           opt->p = opt->p_alloc = 0;
79           opt->h = NULL;
80
81           opt->stopval = -HUGE_VAL;
82           opt->ftol_rel = opt->ftol_abs = 0;
83           opt->xtol_rel = 0; opt->xtol_abs = NULL;
84           opt->maxeval = 0;
85           opt->maxtime = 0;
86           opt->force_stop = 0;
87           opt->force_stop_child = NULL;
88
89           opt->local_opt = NULL;
90           opt->stochastic_population = 0;
91           opt->vector_storage = 0;
92           opt->dx = NULL;
93           opt->work = NULL;
94
95           if (n > 0) {
96                opt->lb = (double *) malloc(sizeof(double) * (n));
97                if (!opt->lb) goto oom;
98                opt->ub = (double *) malloc(sizeof(double) * (n));
99                if (!opt->ub) goto oom;
100                opt->xtol_abs = (double *) malloc(sizeof(double) * (n));
101                if (!opt->xtol_abs) goto oom;
102                nlopt_set_lower_bounds1(opt, -HUGE_VAL);
103                nlopt_set_upper_bounds1(opt, +HUGE_VAL);
104                nlopt_set_xtol_abs1(opt, 0.0);
105           }
106      }
107
108      return opt;
109
110 oom:
111      nlopt_destroy(opt);
112      return NULL;
113 }
114
115 nlopt_opt NLOPT_STDCALL nlopt_copy(const nlopt_opt opt)
116 {
117      nlopt_opt nopt = NULL;
118      unsigned i;
119      if (opt) {
120           nlopt_munge munge;
121           nopt = (nlopt_opt) malloc(sizeof(struct nlopt_opt_s));
122           *nopt = *opt;
123           nopt->lb = nopt->ub = nopt->xtol_abs = NULL;
124           nopt->fc = nopt->h = NULL;
125           nopt->m_alloc = nopt->p_alloc = 0;
126           nopt->local_opt = NULL;
127           nopt->dx = NULL;
128           nopt->work = NULL;
129           opt->force_stop_child = NULL;
130
131           munge = nopt->munge_on_copy;
132           if (munge && nopt->f_data)
133                if (!(nopt->f_data = munge(nopt->f_data))) goto oom;
134
135           if (opt->n > 0) {
136                nopt->lb = (double *) malloc(sizeof(double) * (opt->n));
137                if (!opt->lb) goto oom;
138                nopt->ub = (double *) malloc(sizeof(double) * (opt->n));
139                if (!opt->ub) goto oom;
140                nopt->xtol_abs = (double *) malloc(sizeof(double) * (opt->n));
141                if (!opt->xtol_abs) goto oom;
142                
143                memcpy(nopt->lb, opt->lb, sizeof(double) * (opt->n));
144                memcpy(nopt->ub, opt->ub, sizeof(double) * (opt->n));
145                memcpy(nopt->xtol_abs, opt->xtol_abs, sizeof(double) * (opt->n));
146           }
147
148           if (opt->m) {
149                nopt->m_alloc = opt->m;
150                nopt->fc = (nlopt_constraint *) malloc(sizeof(nlopt_constraint)
151                                                       * (opt->m));
152                if (!nopt->fc) goto oom;
153                memcpy(nopt->fc, opt->fc, sizeof(nlopt_constraint) * (opt->m));
154                for (i = 0; i < opt->m; ++i) nopt->fc[i].tol = NULL;
155                if (munge)
156                     for (i = 0; i < opt->m; ++i)
157                          if (nopt->fc[i].f_data &&
158                              !(nopt->fc[i].f_data
159                                = munge(nopt->fc[i].f_data)))
160                               goto oom;
161                for (i = 0; i < opt->m; ++i)
162                     if (opt->fc[i].tol) {
163                          nopt->fc[i].tol = (double *) malloc(sizeof(double)
164                                                              * nopt->fc[i].m);
165                          if (!nopt->fc[i].tol) goto oom;
166                          memcpy(nopt->fc[i].tol, opt->fc[i].tol,
167                                 sizeof(double) * nopt->fc[i].m);
168                     }
169           }
170
171           if (opt->p) {
172                nopt->p_alloc = opt->p;
173                nopt->h = (nlopt_constraint *) malloc(sizeof(nlopt_constraint)
174                                                      * (opt->p));
175                if (!nopt->h) goto oom;
176                memcpy(nopt->h, opt->h, sizeof(nlopt_constraint) * (opt->p));
177                for (i = 0; i < opt->p; ++i) nopt->h[i].tol = NULL;
178                if (munge)
179                     for (i = 0; i < opt->p; ++i)
180                          if (nopt->h[i].f_data &&
181                              !(nopt->h[i].f_data
182                                = munge(nopt->h[i].f_data)))
183                               goto oom;
184                for (i = 0; i < opt->p; ++i)
185                     if (opt->h[i].tol) {
186                          nopt->h[i].tol = (double *) malloc(sizeof(double)
187                                                              * nopt->h[i].m);
188                          if (!nopt->h[i].tol) goto oom;
189                          memcpy(nopt->h[i].tol, opt->h[i].tol,
190                                 sizeof(double) * nopt->h[i].m);
191                     }
192           }
193
194           if (opt->local_opt) {
195                nopt->local_opt = nlopt_copy(opt->local_opt);
196                if (!nopt->local_opt) goto oom;
197           }
198
199           if (opt->dx) {
200                nopt->dx = (double *) malloc(sizeof(double) * (opt->n));
201                if (!nopt->dx) goto oom;
202                memcpy(nopt->dx, opt->dx, sizeof(double) * (opt->n));
203           }
204      }
205      return nopt;
206
207 oom:
208      nopt->munge_on_destroy = NULL; // better to leak mem than crash
209      nlopt_destroy(nopt);
210      return NULL;
211 }
212
213 /*************************************************************************/
214
215 nlopt_result NLOPT_STDCALL nlopt_set_min_objective(nlopt_opt opt,
216                                                    nlopt_func f, void *f_data)
217 {
218      if (opt) {
219           if (opt->munge_on_destroy) opt->munge_on_destroy(opt->f_data);
220           opt->f = f; opt->f_data = f_data;
221           opt->maximize = 0;
222           if (nlopt_isinf(opt->stopval) && opt->stopval > 0)
223                opt->stopval = -HUGE_VAL; /* switch default from max to min */
224           return NLOPT_SUCCESS;
225      }
226      return NLOPT_INVALID_ARGS;
227 }
228
229 nlopt_result NLOPT_STDCALL nlopt_set_max_objective(nlopt_opt opt, 
230                                                    nlopt_func f, void *f_data)
231 {
232      if (opt) {
233           if (opt->munge_on_destroy) opt->munge_on_destroy(opt->f_data);
234           opt->f = f; opt->f_data = f_data;
235           opt->maximize = 1;
236           if (nlopt_isinf(opt->stopval) && opt->stopval < 0)
237                opt->stopval = +HUGE_VAL; /* switch default from min to max */
238           return NLOPT_SUCCESS;
239      }
240      return NLOPT_INVALID_ARGS;
241 }
242
243 /*************************************************************************/
244
245 nlopt_result
246 NLOPT_STDCALL nlopt_set_lower_bounds(nlopt_opt opt, const double *lb)
247 {
248      if (opt && (opt->n == 0 || lb)) {
249           memcpy(opt->lb, lb, sizeof(double) * (opt->n));
250           return NLOPT_SUCCESS;
251      }
252      return NLOPT_INVALID_ARGS;
253 }
254
255 nlopt_result
256 NLOPT_STDCALL nlopt_set_lower_bounds1(nlopt_opt opt, double lb)
257 {
258      if (opt) {
259           unsigned i;
260           for (i = 0; i < opt->n; ++i)
261                opt->lb[i] = lb;
262           return NLOPT_SUCCESS;
263      }
264      return NLOPT_INVALID_ARGS;
265 }
266
267 nlopt_result
268 NLOPT_STDCALL nlopt_get_lower_bounds(const nlopt_opt opt, double *lb)
269 {
270      if (opt && (opt->n == 0 || lb)) {
271           memcpy(lb, opt->lb, sizeof(double) * (opt->n));
272           return NLOPT_SUCCESS;
273      }
274      return NLOPT_INVALID_ARGS;
275 }
276
277 nlopt_result
278 NLOPT_STDCALL nlopt_set_upper_bounds(nlopt_opt opt, const double *ub)
279 {
280      if (opt && (opt->n == 0 || ub)) {
281           memcpy(opt->ub, ub, sizeof(double) * (opt->n));
282           return NLOPT_SUCCESS;
283      }
284      return NLOPT_INVALID_ARGS;
285 }
286
287 nlopt_result
288 NLOPT_STDCALL nlopt_set_upper_bounds1(nlopt_opt opt, double ub)
289 {
290      if (opt) {
291           unsigned i;
292           for (i = 0; i < opt->n; ++i)
293                opt->ub[i] = ub;
294           return NLOPT_SUCCESS;
295      }
296      return NLOPT_INVALID_ARGS;
297 }
298
299 nlopt_result
300 NLOPT_STDCALL nlopt_get_upper_bounds(const nlopt_opt opt, double *ub)
301 {
302      if (opt && (opt->n == 0 || ub)) {
303           memcpy(ub, opt->ub, sizeof(double) * (opt->n));
304           return NLOPT_SUCCESS;
305      }
306      return NLOPT_INVALID_ARGS;
307 }
308
309 /*************************************************************************/
310
311 #define AUGLAG_ALG(a) ((a) == NLOPT_AUGLAG ||           \
312                        (a) == NLOPT_AUGLAG_EQ ||        \
313                        (a) == NLOPT_LN_AUGLAG ||        \
314                        (a) == NLOPT_LN_AUGLAG_EQ ||     \
315                        (a) == NLOPT_LD_AUGLAG ||        \
316                        (a) == NLOPT_LD_AUGLAG_EQ)
317
318 nlopt_result
319 NLOPT_STDCALL nlopt_remove_inequality_constraints(nlopt_opt opt)
320 {
321      unsigned i;
322      if (!opt) return NLOPT_INVALID_ARGS;
323      if (opt->munge_on_destroy) {
324           nlopt_munge munge = opt->munge_on_destroy;
325           for (i = 0; i < opt->m; ++i)
326                munge(opt->fc[i].f_data);
327      }
328      for (i = 0; i < opt->m; ++i)
329           free(opt->fc[i].tol);
330      free(opt->fc);
331      opt->fc = NULL;
332      opt->m = opt->m_alloc = 0;
333      return NLOPT_SUCCESS;
334 }
335
336 static nlopt_result add_constraint(unsigned *m, unsigned *m_alloc,
337                                    nlopt_constraint **c,
338                                    unsigned fm, nlopt_func fc, nlopt_mfunc mfc,
339                                    void *fc_data,
340                                    const double *tol)
341 {
342      double *tolcopy;
343      unsigned i;
344
345      if ((fc && mfc) || (fc && fm != 1) || (!fc && !mfc))
346           return NLOPT_INVALID_ARGS;
347      if (tol) 
348           for (i = 0; i < fm; ++i) if (tol[i] < 0) return NLOPT_INVALID_ARGS;
349      
350      tolcopy = (double *) malloc(sizeof(double) * fm);
351      if (fm && !tolcopy) return NLOPT_OUT_OF_MEMORY;
352      if (tol)
353           memcpy(tolcopy, tol, sizeof(double) * fm);
354      else
355           for (i = 0; i < fm; ++i) tolcopy[i] = 0;
356
357      *m += 1;
358      if (*m > *m_alloc) {
359           /* allocate by repeated doubling so that 
360              we end up with O(log m) mallocs rather than O(m). */
361           *m_alloc = 2 * (*m);
362           *c = (nlopt_constraint *) realloc(*c,
363                                             sizeof(nlopt_constraint)
364                                             * (*m_alloc));
365           if (!*c) {
366                *m_alloc = *m = 0;
367                free(tolcopy);
368                return NLOPT_OUT_OF_MEMORY;
369           }
370      }
371      
372      (*c)[*m - 1].m = fm;
373      (*c)[*m - 1].f = fc;
374      (*c)[*m - 1].mf = mfc;
375      (*c)[*m - 1].f_data = fc_data;
376      (*c)[*m - 1].tol = tolcopy;
377      return NLOPT_SUCCESS;
378 }
379
380 static int inequality_ok(nlopt_algorithm algorithm) {
381      /* nonlinear constraints are only supported with some algorithms */
382      return (algorithm == NLOPT_LD_MMA || algorithm == NLOPT_LD_CCSAQ 
383              || algorithm == NLOPT_LD_SLSQP
384              || algorithm == NLOPT_LN_COBYLA
385              || AUGLAG_ALG(algorithm) 
386              || algorithm == NLOPT_GN_ISRES
387              || algorithm == NLOPT_GN_ORIG_DIRECT
388              || algorithm == NLOPT_GN_ORIG_DIRECT_L);
389 }
390
391 nlopt_result
392 NLOPT_STDCALL nlopt_add_inequality_mconstraint(nlopt_opt opt, unsigned m,
393                                                nlopt_mfunc fc, void *fc_data,
394                                                const double *tol)
395 {
396      nlopt_result ret;
397      if (!m) { /* empty constraints are always ok */
398           if (opt && opt->munge_on_destroy) opt->munge_on_destroy(fc_data);
399           return NLOPT_SUCCESS;
400      }
401      if (!opt || !inequality_ok(opt->algorithm)) ret = NLOPT_INVALID_ARGS;
402      else ret = add_constraint(&opt->m, &opt->m_alloc, &opt->fc,
403                                m, NULL, fc, fc_data, tol);
404      if (ret < 0 && opt && opt->munge_on_destroy)
405           opt->munge_on_destroy(fc_data);
406      return ret;
407 }
408
409 nlopt_result
410 NLOPT_STDCALL nlopt_add_inequality_constraint(nlopt_opt opt,
411                                               nlopt_func fc, void *fc_data,
412                                               double tol)
413 {
414      nlopt_result ret;
415      if (!opt || !inequality_ok(opt->algorithm)) ret = NLOPT_INVALID_ARGS;
416      else ret = add_constraint(&opt->m, &opt->m_alloc, &opt->fc,
417                                1, fc, NULL, fc_data, &tol);
418      if (ret < 0 && opt && opt->munge_on_destroy)
419           opt->munge_on_destroy(fc_data);
420      return ret;
421 }
422
423 nlopt_result
424 NLOPT_STDCALL nlopt_remove_equality_constraints(nlopt_opt opt)
425 {
426      unsigned i;
427      if (!opt) return NLOPT_INVALID_ARGS;
428      if (opt->munge_on_destroy) {
429           nlopt_munge munge = opt->munge_on_destroy;
430           for (i = 0; i < opt->p; ++i)
431                munge(opt->h[i].f_data);
432      }
433      for (i = 0; i < opt->p; ++i)
434           free(opt->h[i].tol);
435      free(opt->h);
436      opt->h = NULL;
437      opt->p = opt->p_alloc = 0;
438      return NLOPT_SUCCESS;
439 }
440
441 static int equality_ok(nlopt_algorithm algorithm) {
442      /* equality constraints (h(x) = 0) only via some algorithms */
443      return (AUGLAG_ALG(algorithm) 
444              || algorithm == NLOPT_LD_SLSQP
445              || algorithm == NLOPT_GN_ISRES
446              || algorithm == NLOPT_LN_COBYLA);
447 }
448
449 nlopt_result
450 NLOPT_STDCALL nlopt_add_equality_mconstraint(nlopt_opt opt, unsigned m,
451                                              nlopt_mfunc fc, void *fc_data,
452                                              const double *tol)
453 {
454      nlopt_result ret;
455      if (!m) { /* empty constraints are always ok */
456           if (opt && opt->munge_on_destroy) opt->munge_on_destroy(fc_data);
457           return NLOPT_SUCCESS;
458      }
459      if (!opt || !equality_ok(opt->algorithm)
460          || nlopt_count_constraints(opt->p, opt->h) + m > opt->n) 
461           ret = NLOPT_INVALID_ARGS;
462      else ret = add_constraint(&opt->p, &opt->p_alloc, &opt->h,
463                                m, NULL, fc, fc_data, tol);
464      if (ret < 0 && opt && opt->munge_on_destroy)
465           opt->munge_on_destroy(fc_data);
466      return ret;
467 }
468
469 nlopt_result
470 NLOPT_STDCALL nlopt_add_equality_constraint(nlopt_opt opt,
471                                             nlopt_func fc, void *fc_data,
472                                             double tol)
473 {
474      nlopt_result ret;
475      if (!opt || !equality_ok(opt->algorithm)
476          || nlopt_count_constraints(opt->p, opt->h) + 1 > opt->n)
477           ret = NLOPT_INVALID_ARGS;
478      else ret = add_constraint(&opt->p, &opt->p_alloc, &opt->h,
479                                1, fc, NULL, fc_data, &tol);
480      if (ret < 0 && opt && opt->munge_on_destroy)
481           opt->munge_on_destroy(fc_data);
482      return ret;
483 }
484
485 /*************************************************************************/
486
487 #define SET(param, T, arg)                                              \
488    nlopt_result NLOPT_STDCALL nlopt_set_##param(nlopt_opt opt, T arg)   \
489    {                                                                    \
490         if (opt) {                                                      \
491              opt->arg = arg;                                            \
492              return NLOPT_SUCCESS;                                      \
493         }                                                               \
494         return NLOPT_INVALID_ARGS;                                      \
495    }
496
497
498 #define GET(param, T, arg) T NLOPT_STDCALL      \
499    nlopt_get_##param(const nlopt_opt opt) {     \
500         return opt->arg;                        \
501    }
502
503 #define GETSET(param, T, arg) GET(param, T, arg) SET(param, T, arg)
504
505 GETSET(stopval, double, stopval)
506
507 GETSET(ftol_rel, double, ftol_rel)
508 GETSET(ftol_abs, double, ftol_abs)
509 GETSET(xtol_rel, double, xtol_rel)
510
511 nlopt_result
512 NLOPT_STDCALL nlopt_set_xtol_abs(nlopt_opt opt, const double *xtol_abs)
513 {
514      if (opt) {
515           memcpy(opt->xtol_abs, xtol_abs, opt->n * sizeof(double));
516           return NLOPT_SUCCESS;
517      }
518      return NLOPT_INVALID_ARGS;
519 }
520
521 nlopt_result
522 NLOPT_STDCALL nlopt_set_xtol_abs1(nlopt_opt opt, double xtol_abs)
523 {
524      if (opt) {
525           unsigned i;
526           for (i = 0; i < opt->n; ++i)
527                opt->xtol_abs[i] = xtol_abs;
528           return NLOPT_SUCCESS;
529      }
530      return NLOPT_INVALID_ARGS;
531 }
532
533 nlopt_result
534 NLOPT_STDCALL nlopt_get_xtol_abs(const nlopt_opt opt, double *xtol_abs)
535 {
536      memcpy(xtol_abs, opt->xtol_abs, opt->n * sizeof(double));
537      return NLOPT_SUCCESS;
538 }
539
540 GETSET(maxeval, int, maxeval)
541
542 GETSET(maxtime, double, maxtime)
543
544 /*************************************************************************/
545
546 nlopt_result
547 NLOPT_STDCALL nlopt_set_force_stop(nlopt_opt opt, int force_stop)
548 {
549      if (opt) {
550           opt->force_stop = force_stop;
551           if (opt->force_stop_child)
552                return nlopt_set_force_stop(opt->force_stop_child, force_stop);
553           return NLOPT_SUCCESS;
554      }
555      return NLOPT_INVALID_ARGS;
556 }
557
558 GET(force_stop, int, force_stop)
559 nlopt_result NLOPT_STDCALL nlopt_force_stop(nlopt_opt opt) { 
560      return nlopt_set_force_stop(opt, 1); 
561 }
562
563 /*************************************************************************/
564
565 GET(algorithm, nlopt_algorithm, algorithm)
566 GET(dimension, unsigned, n)
567
568 /*************************************************************************/
569
570 nlopt_result
571 NLOPT_STDCALL nlopt_set_local_optimizer(nlopt_opt opt,
572                                         const nlopt_opt local_opt)
573 {
574      if (opt) {
575           if (local_opt && local_opt->n != opt->n) return NLOPT_INVALID_ARGS;
576           nlopt_destroy(opt->local_opt);
577           opt->local_opt = nlopt_copy(local_opt);
578           if (local_opt) {
579                if (!opt->local_opt) return NLOPT_OUT_OF_MEMORY;
580                nlopt_set_lower_bounds(opt->local_opt, opt->lb);
581                nlopt_set_upper_bounds(opt->local_opt, opt->ub);
582                nlopt_remove_inequality_constraints(opt->local_opt);
583                nlopt_remove_equality_constraints(opt->local_opt);
584                nlopt_set_min_objective(opt->local_opt, NULL, NULL);
585                nlopt_set_munge(opt->local_opt, NULL, NULL);
586                opt->local_opt->force_stop = 0;
587           }
588           return NLOPT_SUCCESS;
589      }
590      return NLOPT_INVALID_ARGS;
591 }
592
593 /*************************************************************************/
594
595 GETSET(population, unsigned, stochastic_population)
596 GETSET(vector_storage, unsigned, vector_storage)
597
598 /*************************************************************************/
599
600 nlopt_result NLOPT_STDCALL nlopt_set_initial_step1(nlopt_opt opt, double dx)
601 {
602      unsigned i;
603      if (!opt || dx == 0) return NLOPT_INVALID_ARGS;
604      if (!opt->dx && opt->n > 0) {
605           opt->dx = (double *) malloc(sizeof(double) * (opt->n));
606           if (!opt->dx) return NLOPT_OUT_OF_MEMORY;
607      }
608      for (i = 0; i < opt->n; ++i) opt->dx[i] = dx;
609      return NLOPT_SUCCESS;
610 }
611
612 nlopt_result
613 NLOPT_STDCALL nlopt_set_initial_step(nlopt_opt opt, const double *dx)
614 {
615      unsigned i;
616      if (!opt) return NLOPT_INVALID_ARGS;
617      if (!dx) {
618           free(opt->dx); opt->dx = NULL;
619           return NLOPT_SUCCESS;
620      }
621      for (i = 0; i < opt->n; ++i) if (dx[i] == 0) return NLOPT_INVALID_ARGS;
622      if (!opt->dx && nlopt_set_initial_step1(opt, 1) == NLOPT_OUT_OF_MEMORY)
623           return NLOPT_OUT_OF_MEMORY;
624      memcpy(opt->dx, dx, sizeof(double) * (opt->n));
625      return NLOPT_SUCCESS;
626 }
627
628 nlopt_result
629 NLOPT_STDCALL nlopt_get_initial_step(const nlopt_opt opt, const double *x, 
630                                      double *dx)
631 {
632      if (!opt) return NLOPT_INVALID_ARGS;
633      if (!opt->n) return NLOPT_SUCCESS;
634      if (!opt->dx) {
635           nlopt_opt o = (nlopt_opt) opt; /* discard const temporarily */
636           nlopt_result ret = nlopt_set_default_initial_step(o, x);
637           if (ret != NLOPT_SUCCESS) return ret;
638           memcpy(dx, o->dx, sizeof(double) * (opt->n));
639           free(o->dx); o->dx = NULL; /* don't save, since x-dependent */
640      }
641      else
642           memcpy(dx, opt->dx, sizeof(double) * (opt->n));
643      return NLOPT_SUCCESS;
644 }
645
646 nlopt_result
647 NLOPT_STDCALL nlopt_set_default_initial_step(nlopt_opt opt, const double *x)
648 {
649      const double *lb, *ub;
650      unsigned i;
651
652      if (!opt || !x) return NLOPT_INVALID_ARGS;
653      lb = opt->lb; ub = opt->ub;
654
655      if (!opt->dx && nlopt_set_initial_step1(opt, 1) == NLOPT_OUT_OF_MEMORY)
656           return NLOPT_OUT_OF_MEMORY;
657
658      /* crude heuristics for initial step size of nonderivative algorithms */
659      for (i = 0; i < opt->n; ++i) {
660           double step = HUGE_VAL;
661
662           if (!nlopt_isinf(ub[i]) && !nlopt_isinf(lb[i])
663               && (ub[i] - lb[i]) * 0.25 < step && ub[i] > lb[i])
664                step = (ub[i] - lb[i]) * 0.25;
665           if (!nlopt_isinf(ub[i]) 
666               && ub[i] - x[i] < step && ub[i] > x[i])
667                step = (ub[i] - x[i]) * 0.75;
668           if (!nlopt_isinf(lb[i]) 
669               && x[i] - lb[i] < step && x[i] > lb[i])
670                step = (x[i] - lb[i]) * 0.75;
671
672           if (nlopt_isinf(step)) {
673                if (!nlopt_isinf(ub[i]) 
674                    && fabs(ub[i] - x[i]) < fabs(step))
675                     step = (ub[i] - x[i]) * 1.1;
676                if (!nlopt_isinf(lb[i]) 
677                    && fabs(x[i] - lb[i]) < fabs(step))
678                     step = (x[i] - lb[i]) * 1.1;
679           }
680           if (nlopt_isinf(step) || step == 0) {
681                step = x[i];
682           }
683           if (nlopt_isinf(step) || step == 0)
684                step = 1;
685           
686           opt->dx[i] = step;
687      }
688      return NLOPT_SUCCESS;
689 }
690
691 /*************************************************************************/
692
693 void NLOPT_STDCALL nlopt_set_munge(nlopt_opt opt,
694                                    nlopt_munge munge_on_destroy,
695                                    nlopt_munge munge_on_copy) {
696      if (opt) {
697           opt->munge_on_destroy = munge_on_destroy;
698           opt->munge_on_copy = munge_on_copy;
699      }
700 }
701
702 /*************************************************************************/