TPCCLIB
Loading...
Searching...
No Matches
Functions
praxis.c File Reference

Powell-Brent (Praxis) nonlinear optimization. More...

#include "tpcclibConfig.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
#include "tpcextensions.h"
#include "tpcrand.h"
#include "tpcnlopt.h"

Go to the source code of this file.

Functions

int nloptPowellBrent (NLOPT *nlo, unsigned int ktm, TPCSTATUS *status)
 Powell-Brent (Praxis) non-linear unconstrained optimization.

Detailed Description

Powell-Brent (Praxis) nonlinear optimization.

Note
Not yet for production use.
Todo
Test with practical applications.

Definition in file praxis.c.

Function Documentation

◆ nloptPowellBrent()

int nloptPowellBrent ( NLOPT * nlo,
unsigned int ktm,
TPCSTATUS * status )

Powell-Brent (Praxis) non-linear unconstrained optimization.

Praxis is a general purpose routine for the minimization of a function in several variables. The algorithm used is a modification of conjugate gradient search method by Powell, with changes by Brent, who gave an algol-w program, which served as a basis for this function.

References:

  1. Powell, MJD (1964) An efficient method for finding the minimum of a function in several variables without calculating derivatives. Computer Journal 7:155-162.
  2. Brent, RP (1973) Algorithms for minimization without derivatives. Prentice Hall, Englewood Cliffs. Dover 2013 republication.
Returns
enum tpcerror (TPCERROR_OK when successful).
Author
Vesa Oikonen
Parameters
nloPointer to NLOPT struct. Initial guess must be given in x[]. Initial step sizes must be given in xdelta[]. Note that constraints are not applied here. Parameter tolerances are used as stopping criteria.
ktmParameter tolerances (xtol[]) are given in the previous struct. Praxis stops if the criterion 2 * ||x[k]-x[k-1]|| <= sqrt(macheps) * ||x[k]|| + xtol is fulfilled more than ktm times for every parameter x. You should first try with ktm=1, and usually not higher than 4.
statusPointer to status data; enter NULL if not needed.

Definition at line 59 of file praxis.c.

74 {
75 int verbose=0; if(status!=NULL) verbose=status->verbose;
76 if(verbose>0) printf("%s()\n", __func__);
77 if(nlo==NULL) {
78 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
79 return TPCERROR_FAIL;
80 }
81 if(nlo->totalNr<1 || nlo->xfull==NULL || nlo->_fun==NULL) {
82 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
83 return TPCERROR_NO_DATA;
84 }
85
86
87 /* Check deltas */
88 unsigned int fixedNr=0;
89 for(unsigned int i=0; i<nlo->totalNr; i++) if(nlo->xdelta[i]<=0.0) fixedNr++;
90 if(fixedNr==nlo->totalNr) {
91 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_PARNR);
92 return TPCERROR_INVALID_PARNR;
93 }
94
95
96 /* Set a safe machine epsilon and related values */
97 double macheps=doubleMachEps()*100.0;
98 if(verbose>1) printf("macheps := %g\n", macheps);
99 double m2=sqrt(macheps);
100 double m4=sqrt(m2);
101 if(verbose>1) {
102 printf("m2 := %e\n", m2);
103 printf("m4 := %e\n", m4);
104 }
105 double small=macheps*macheps;
106 double vsmall=small*small;
107 double large=1.0/small;
108 double vlarge=1.0/vsmall;
109 if(verbose>1) {
110 printf("small := %e\n", small);
111 printf("vsmall := %e\n", vsmall);
112 printf("large := %e\n", large);
113 printf("vlarge := %e\n", vlarge);
114 }
115
116 /* Set common initial step size and scaling parameter */
117 double h=0.0, scbd=1.0, t=0.0, t2;
118 {
119 unsigned int i;
120 double min=nan(""), max=nan("");
121 for(i=0; i<nlo->totalNr; i++) if(nlo->xdelta[i]>0.0) {
122 h+=nlo->xdelta[i];
123 t+=nlo->xtol[i];
124 if(isnan(min) || nlo->xdelta[i]<min) min=nlo->xdelta[i];
125 if(isnan(max) || nlo->xdelta[i]>max) max=nlo->xdelta[i];
126 }
127 h/=(double)(nlo->totalNr-fixedNr);
128 t/=(double)(nlo->totalNr-fixedNr);
129 if(verbose>2) {printf("xdelta_min := %g\nxdelta_max := %g\n", min, max); fflush(stdout);}
130 scbd+=log10(max/min); // >1 if param scales are different
131 t2=small+fabs(t); t=t2;
132 if(h<100.0*t) h=100.0*t;
133 }
134 if(verbose>1) {
135 printf("step := %g\n", h);
136 printf("scbd := %g\n", scbd);
137 printf("t(2) := %g\n", t);
138 fflush(stdout);
139 }
140 if(isnan(scbd)) { // checking that user provided xdeltas
141 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
142 return TPCERROR_INVALID_VALUE;
143 }
144
145
146 /* Is the problem ill-conditioned (1) or not?
147 This variable is automatically set, when Praxis finds the problem to be
148 ill-conditioned during iterations. */
149 int illc=0;
150
151 /* Initialize */
152 if(verbose>2) {printf("initializing\n"); fflush(stdout);}
153 unsigned int i, j, dim=nlo->totalNr;
154 double ldfac; if(illc) ldfac=0.1; else ldfac=0.01;
155 unsigned int nl=0, kt=0;
156 unsigned int nf=0; // nr of objective function calls
157 double dmin=small;
158 double ldt=h;
159 double d[dim], y[dim], z[dim];
160 // double v[dim][dim];
161 double **v=malloc(dim*sizeof(double *)); if(v==NULL) {
162 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
163 return TPCERROR_OUT_OF_MEMORY;
164 }
165 for(i=0; i<dim; i++) {
166 v[i]=malloc(dim*sizeof(double));
167 if(v[i]==NULL) {
168 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
169 return TPCERROR_OUT_OF_MEMORY;
170 }
171 }
172 /* Make a copy of parameter list */
173 if(verbose>2) {printf("copying parameters\n"); fflush(stdout);}
174 double p[dim];
175 for(i=0; i<dim; i++) p[i]=nlo->xfull[i];
176
177 for(i=0; i<dim; i++) {
178 d[i]=y[i]=z[i]=0.0;
179 for(j=0; j<dim; j++) v[i][j]=(i==j ? 1.0 : 0.0);
180 }
181 double qd0=0.0, q0[dim], q1[dim];
182 for(i=0; i<dim; i++) q1[i]=q0[i]=nlo->xfull[i];
183
184
185 /* Calculate function value with the initial guess */
186 if(verbose>2) {printf("evaluating initial guess\n"); fflush(stdout);}
187 double fx, qf1;
188 fx=qf1=(*nlo->_fun)(nlo->totalNr, nlo->xfull, nlo->fundata);
189 nf++;
190 if(!isfinite(fx)) {
191 for(i=0; i<dim; i++) free(v[i]);
192 free(v);
193 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
194 return TPCERROR_INVALID_VALUE;
195 }
196 if(verbose>1) {printf("initial_objfunc := %g\n", fx); fflush(stdout);}
197
198
199 /* The main loop */
200 double sf, s, dn, qd1=0.0;
201 int done=0, nLoop=0, ret=0;
202 while(!done) {
203 nLoop++;
204 if(verbose>2) {
205 printf("\n-------------------------------\nloop %d, with %d fcalls\n", nLoop, nf);
206 fflush(stdout);
207 }
208 sf=d[0]; s=d[0]=0.0;
209
210 /* Minimize along the first direction */
211 ret=praxisMin(0, 2, &d[0], &s, fx, 0, &fx, v, dim, h, t,
212 macheps, m2, m4, small, dmin, ldt, &nl,
213 &nf, qd0, qd1, q0, q1,
214 p, nlo->_fun, nlo->fundata);
215 if(ret) {
216 if(verbose>0) printf("Error %d in praxisMin()\n", ret);
217 for(i=0; i<dim; i++) free(v[i]);
218 free(v);
219 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
220 return TPCERROR_FAIL;
221 }
222
223 if((sf<=(0.9*d[0])) || ((0.9*sf)>=d[0])) for(i=1; i<dim; i++) d[i]=0.0;
224
225 /* Other directions */
226 unsigned int k, kl=0, k2;
227 double sl, df, f1, lds;
228 for(k=1; k<dim; k++) {
229 if(verbose>5) printf("direction %d\n", k);
230 for(i=0; i<dim; i++) y[i]=p[i];
231 sf=fx; if(kt>0) illc=1;
232 do { // Usually only once, but twice if noticed here that ill-conditioned
233 kl=k; df=0.0;
234 /* If ill-conditioned, take random step to get off resolution valley */
235 if(illc) {
236 for(i=0; i<dim; i++) {
237 z[i]=(0.1*ldt+t2*pow(10.0,(double)kt)) * (drand()-0.5);
238 s=z[i]; for(j=0; j<dim; j++) p[j]+=s*v[j][i];
239 }
240 /* Calculate function value */
241 fx=(*nlo->_fun)(dim, p, nlo->fundata);
242 if(!isfinite(fx)) {
243 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
244 return TPCERROR_INVALID_VALUE;
245 }
246 nf++;
247 }
248 /* Minimize along non-conjugate directions */
249 for(k2=k; k2<dim; k2++) {
250 if(verbose>6) printf(" non-conjugate direction %d\n", k2);
251 sl=fx; s=0.0;
252 ret=praxisMin(k2, 2, &d[k2], &s, fx, 0, &fx, v, dim, h, t,
253 macheps, m2, m4, small, dmin, ldt, &nl,
254 &nf, qd0, qd1, q0, q1,
255 p, nlo->_fun, nlo->fundata);
256 if(ret) {
257 if(verbose>0) printf("Error %d in praxisMin()\n", ret);
258 for(i=0; i<dim; i++) free(v[i]);
259 free(v);
260 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
261 return TPCERROR_FAIL;
262 }
263 if(illc) s=d[k2]*(s+z[k2])*(s+z[k2]); else s=sl-fx;
264 if(df<s) {df=s; kl=k2;}
265 }
266 if(illc) break;
267 if(df>=fabs(100.0*macheps*fx)) break;
268 illc=1;
269 } while(1);
270
271 /* Minimize along conjugate directions */
272 for(k2=0; k2<k; k2++) { // WHY NOT UNTIL k2==k ??????
273// for(k2=0; k2<=k; k2++) {
274 if(verbose>6) printf(" conjugate direction %d\n", k2);
275 s=0.0;
276 ret=praxisMin(k2, 2, &d[k2], &s, fx, 0, &fx, v, dim, h, t,
277 macheps, m2, m4, small, dmin, ldt, &nl,
278 &nf, qd0, qd1, q0, q1,
279 p, nlo->_fun, nlo->fundata);
280 if(ret) {
281 if(verbose>0) printf("Error %d in praxisMin()\n", ret);
282 for(i=0; i<dim; i++) free(v[i]);
283 free(v);
284 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
285 return TPCERROR_FAIL;
286 }
287 }
288 f1=fx; fx=sf; lds=0.0;
289 for(i=0; i<dim; i++) {
290 sl=p[i]; p[i]=y[i]; y[i]=sl-y[i]; sl=y[i]; lds+=sl*sl;
291 }
292 lds=sqrt(lds);
293 if(lds>small) {
294 for(i=kl-1; i>=k; i--) {
295 for(j=0; j<dim; j++) v[j][i+1]=v[j][i];
296 d[i+1]=d[i];
297 }
298 d[k]=0.0; for(i=0; i<dim; i++) v[i][k]=y[i]/lds;
299 ret=praxisMin(k, 4, &d[k], &lds, f1, 1, &fx, v, dim, h, t,
300 macheps, m2, m4, small, dmin, ldt, &nl,
301 &nf, qd0, qd1, q0, q1,
302 p, nlo->_fun, nlo->fundata);
303 if(ret) {
304 if(verbose>0) printf("Error %d in praxisMin()\n", ret);
305 for(i=0; i<dim; i++) free(v[i]);
306 free(v);
307 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
308 return TPCERROR_FAIL;
309 }
310 if(lds<=0.0) {lds=-lds; for(i=0; i<dim; i++) v[i][k]=-v[i][k];}
311 }
312 ldt*=ldfac; if(ldt<lds) ldt=lds;
313 for(i=0, t2=0.0; i<dim; i++) t2+=p[i]*p[i];
314 t2=sqrt(t2); t2*=m2; t2+=t;
315
316 if(verbose>1) printf(" objfunc := %g\n", fx);
317
318 /* Stopping criterion */
319 if(ldt>(0.5*t2)) kt=0; else kt++;
320 if(verbose>3 && kt>0) printf(" kt := %d\n", kt);
321 if(kt>ktm) {done=1; break;}
322
323 } // next direction k
324 if(done) break;
325
326 /*
327 * We are probably in a curved valley, try quadratic extrapolation
328 */
329 /* but only if we have more than one parameter to fit */
330 if(dim<2) {
331 if(nf>1000) done=1;
332 continue;
333 }
334
335 /* Praxis quad:
336 looks for the minimum along a curve defined by q0, q1, and x (q2). */
337 {
338 double qa, qb, qc, qs, ql;
339 qs=fx; fx=qf1; qf1=qs;
340 for(i=0, qd1=0.0; i<dim; i++) {
341 qs=p[i]; ql=q1[i]; p[i]=ql; q1[i]=qs;
342 qd1+=(qs-ql)*(qs-ql);
343 }
344 qd1=sqrt(qd1); ql=qd1; qs=0.0;
345 if(qd0>0.0 && qd1>0.0 && nl>=3*dim*dim) {
346 ret=praxisMin(-1, 2, &qs, &ql, qf1, 1, &fx, v, dim, h, t,
347 macheps, m2, m4, small, dmin, ldt, &nl,
348 &nf, qd0, qd1, q0, q1,
349 p, nlo->_fun, nlo->fundata);
350 if(ret) {
351 if(verbose>0) printf("Error %d in praxisMin()\n", ret);
352 for(i=0; i<dim; i++) free(v[i]);
353 free(v);
354 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
355 return TPCERROR_FAIL;
356 }
357 qa = ql*(ql-qd1)/(qd0*(qd0+qd1));
358 qb = (ql+qd0)*(qd1-ql)/(qd0*qd1);
359 qc = ql*(ql+qd0)/(qd1*(qd0+qd1));
360 } else {
361 fx=qf1; qa=qb=0.0; qc=1.0;
362 }
363 if(!isfinite(qa)) qa=0.0;
364 if(!isfinite(qb)) qb=0.0;
365 if(!isfinite(qc)) qc=1.0;
366 qd0=qd1;
367 for(i=0; i<dim; i++) {
368 qs=q0[i]; q0[i]=p[i];
369 p[i]=qa*s + qb*p[i] + qc*q1[i];
370 }
371 }
372 dn=0.0;
373 for(i=0; i<dim; i++) {d[i]=1.0/sqrt(d[i]); if(dn<d[i]) dn=d[i];}
374 for(j=0; j<dim; j++) {
375 s=d[j]/dn;
376 if(isfinite(s)) for(i=0; i<dim; i++) v[i][j]*=s;
377 }
378 /* Scale axis to reduce condition number */
379 if(scbd>1.0) {
380 s=vlarge;
381 for(i=0; i<dim; i++) {
382 sl=0.0;
383 for(j=0; j<dim; j++) sl+=v[i][j]*v[i][j];
384 z[i]=sqrt(sl); if(!isfinite(z[i]) || z[i]<m4) z[i]=m4;
385 if(s>z[i]) s=z[i];
386 }
387 for(i=0; i<dim; i++) {
388 sl=s/z[i]; z[i]=1.0/sl;
389 if(z[i]>scbd) {sl=1.0/scbd; z[i]=scbd;}
390 for(j=0; j<dim; j++) v[i][j]*=sl;
391 }
392 }
393 /* Calculate new set of orthogonal directions before repeating
394 the main loop; first, transpose v[][] for the Praxis minfit */
395 for(i=1; i<dim; i++) for(j=0; j<i; j++) {
396 s=v[i][j]; v[i][j]=v[j][i]; v[j][i]=s;
397 }
398 /* Call Praxis minfit to find the SVD of v[][].
399 This gives the principal values and principal directions of approximating
400 quadratic form without squaring the condition number */
401 ret=praxisMinfit((int)dim, macheps, vsmall, v, d);
402 if(ret) {
403 if(verbose>0) printf("Error %d in praxisMinfit()\n", ret);
404 for(i=0; i<dim; i++) free(v[i]);
405 free(v);
406 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
407 return TPCERROR_FAIL;
408 }
409 if(scbd>1.0) {
410 for(i=0; i<dim; i++) {s=z[i]; for(j=0; j<dim; j++) v[i][j]*=s;}
411 for(i=0; i<dim; i++) {
412 s=0.0; for(j=0; j<dim; j++) s+=v[j][i]*v[j][i];
413 s=sqrt(s); d[i]*=s; s=1.0/s; for(j=0; j<dim; j++) v[j][i]*=s;
414 }
415 }
416 for(i=0; i<dim; i++) {
417 if((dn*d[i]) > large) d[i]=vsmall; // dn was calculated above
418 else if((dn*d[i]) < small) d[i]=vlarge;
419 else d[i]=pow(dn*d[i], -2.0);
420 }
421 /* the new eigenvalues and eigenvectors */
422 praxisSort(d, v, dim);
423 dmin=d[dim-1]; if(dmin<small) dmin=small;
424 illc=(m2*d[0]) > dmin;
425
426
427
428
429 if(nf>1000) done=1;
430 }
431 if(verbose>1) printf("func_evals := %d\n", nf);
432
433
434 /* Free allocated memory */
435 for(i=0; i<dim; i++) free(v[i]);
436 free(v);
437
438 /* If ok, then copy parameters to the provided struct */
439 for(i=0; i<dim; i++) nlo->xfull[i]=p[i];
440
441
442 // Just to prevent warning for now
443 if(sf>0.0) sf=0.0;
444 if(kt>0) kt=0;
445 if(nl>0) nl=0;
446 if(ldfac>0.0) ldfac=0.0;
447 if(ktm>0) ktm=0;
448
449
450 return TPCERROR_OK;
451}
doubleMachEps
double doubleMachEps()
Definition doubleutil.c:105
drand
double drand()
Definition gaussdev.c:66
statusSet
void statusSet(TPCSTATUS *s, const char *func, const char *srcfile, int srcline, tpcerror error)
Definition statusmsg.c:142
NLOPT::_fun
double(* _fun)(int, double *, void *)
Definition tpcnlopt.h:42
NLOPT::fundata
void * fundata
Definition tpcnlopt.h:44
NLOPT::xfull
double * xfull
Definition tpcnlopt.h:29
NLOPT::xdelta
double * xdelta
Definition tpcnlopt.h:36
NLOPT::xtol
double * xtol
Definition tpcnlopt.h:39
NLOPT::totalNr
unsigned int totalNr
Definition tpcnlopt.h:27
TPCSTATUS::verbose
int verbose
Verbose level, used by statusPrint() etc.
Definition tpcextensions.h:242
TPCERROR_INVALID_VALUE
@ TPCERROR_INVALID_VALUE
Invalid value.
Definition tpcextensions.h:220
TPCERROR_INVALID_PARNR
@ TPCERROR_INVALID_PARNR
Invalid number of parameters.
Definition tpcextensions.h:222
TPCERROR_FAIL
@ TPCERROR_FAIL
General error.
Definition tpcextensions.h:192
TPCERROR_OUT_OF_MEMORY
@ TPCERROR_OUT_OF_MEMORY
Cannot allocate memory.
Definition tpcextensions.h:193
TPCERROR_OK
@ TPCERROR_OK
No error.
Definition tpcextensions.h:191
TPCERROR_NO_DATA
@ TPCERROR_NO_DATA
File contains no data.
Definition tpcextensions.h:203
Generated on for TPCCLIB by doxygen 1.15.0