powell.c File Reference

Powell function minimization routines. More...

#include "libtpcmodel.h"

Go to the source code of this file.

Functions
int	powell (double *p, double *delta, int parNr, double ftol, int *iterNr, double *fret, double(*_fun)(int, double *, void *), void *fundata, int verbose)

Variables
int	POWELL_LINMIN_MAXIT =100

Detailed Description

Powell function minimization routines.

Author

Vesa Oikonen

Based on Numerical recipes in C (Press et al.).

Definition in file powell.c.

Function Documentation

powell()

int powell	(	double *	p,
		double *	delta,
		int	parNr,
		double	ftol,
		int *	iterNr,
		double *	fret,
		double(*	_fun )(int, double *, void *),
		void *	fundata,
		int	verbose )

Powell function minimization routine.

Returns

Returns 0, if successful, 1 if required tolerance was not reached, 2 if initial guess does not give finite function value, 3 if final function value is NaN or infinite, and >3 in case of another error.

See also

simplex, tgo, bobyqa, nlopt1D

Parameters

p	Initial guess and final set of parameters
delta	Initial changes for parameters, ==0 if fixed
parNr	Nr of parameters
ftol	Fractional tolerance (for WSS); 0<ftol<1
iterNr	Max nr of iterations, and nr of required iters
fret	Function return value (WSS) at minimum
_fun	Function to minimize (must return the WSS)
fundata	Pointer to data which is passed on to the function; NULL if not needed
verbose	Verbose level; if zero, then nothing is printed into stdout or stderr

Definition at line 43 of file powell.c.

  {
  int pbIterNr;
  int i, j, ibig, iterMax, fixed[MAX_PARAMETERS];
  double xi[MAX_PARAMETERS][MAX_PARAMETERS]; /* Matrix for directions */
  double pt[MAX_PARAMETERS], ptt[MAX_PARAMETERS], xit[MAX_PARAMETERS];
  double del, fp, fptt, t;
  double origp[MAX_PARAMETERS];
  int ftol_reached=0;
 
 
  if(verbose>0) printf("in powell(,,%d,%g,%d,,)\n", parNr, ftol, *iterNr);
  if(verbose>1) {
    printf("Initial parameter guesses and deltas:\n");
    for(i=0; i<parNr; i++) printf("  %g  %g\n", p[i], delta[i]);
  }
  *fret=nan("");
  if(p==NULL) return(11);
  if(delta==NULL) return(12);
  if(parNr<1) return(21);
  if(ftol<=0.0) return(22);
  if(ftol>=1.0) return(23);
  if((*iterNr)<1) return(24);
 
  /* SetUp */
  _powellFunc=_fun;
  _powellFuncData=fundata;
  iterMax=*iterNr; /* save the max nr of iterations */
  _powell_ncom=parNr;
  /* Function value at initial point */
  _powell_func_calls=1;
  *fret=(*_powellFunc)(parNr, p, _powellFuncData);
  if(verbose>10) printf("initial point fret=%g\n", *fret);
  if(!isfinite(*fret)) {
    if(verbose>0) printf("in powell(): objf failed at initial point.\n");
    *fret=nan(""); return(2);
  }
  /* Save the initial point (pt[] will be changed later) */
  for(j=0; j<parNr; j++) origp[j]=pt[j]=p[j];
  /* Check which parameters are fixed */
  for(i=0; i<parNr; i++) if(fabs(delta[i])<1.0e-20) fixed[i]=1; else fixed[i]=0;
 
  /* Initiate matrix for directions */
  for(i=0; i<parNr; i++)
    for(j=0; j<parNr; j++)
      if(i==j) xi[i][j]=delta[i]; else xi[i][j]=0.0;
 
  /* Iterate */
  for(*iterNr=1; ; (*iterNr)++) {
    if(verbose>2) printf("  iteration %d\n", *iterNr);
    fp=*fret; ibig=0; del=0.0; /* largest function decrease */
 
    /* In each iteration, loop over all directions in the set */
    for(i=0; i<parNr; i++) {
      if(fixed[i]) continue; /* do nothing with fixed parameters */
      for(j=0; j<parNr; j++) if(fixed[j]) xit[j]=0.0; else xit[j]=xi[j][i];
      fptt=*fret;
      /* minimize along direction xit */
      pbIterNr=POWELL_LINMIN_MAXIT;
      _powell_linmin(p, xit, parNr, fret, &pbIterNr);
      if(verbose>3) printf("iterNr in _powell_linmin() with p%d: %d\n",
                           i, pbIterNr);
      if(fabs(fptt-(*fret))>del) {del=fabs(fptt-(*fret)); ibig=i;}
    }
    if(verbose>20) printf("fret=%g  fp=%g\n", *fret, fp);
 
    /* Check if done */
#if(0)
    if(2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) break;
#else
    if(2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) {
      if(ftol_reached>0 || (*iterNr)>=iterMax) break; else ftol_reached++;
    } else ftol_reached=0;
#endif
    if((*iterNr)>=iterMax) {
      if(verbose>0) printf("max iterations nr exceeded in powell().\n");
      break;
    }
    /* Construct the extrapolated point and the average direction moved */
    for(j=0; j<parNr; j++) {
      ptt[j]=2.0*p[j]-pt[j]; xit[j]=p[j]-pt[j];
      pt[j]=p[j]; /* save the old starting point */
    }
    fptt=(*_powellFunc)(parNr, ptt, _powellFuncData); _powell_func_calls++;
    if(fptt<fp) {
      t=2.0*(fp-2.0*(*fret)+fptt)*_powell_sqr(fp-(*fret)-del)-del*_powell_sqr(fp-fptt);
      if(t<0.0) {
        pbIterNr=POWELL_LINMIN_MAXIT;
        _powell_linmin(p, xit, parNr, fret, &pbIterNr);
        if(verbose>3) printf("iterNr in _powell_linmin(): %d\n", pbIterNr);
        for(j=0; j<parNr; j++) {
          xi[j][ibig]=xi[j][parNr-1]; xi[j][parNr-1]=xit[j];}
      }
    }
  } /* next iteration */
  if(verbose>1) {
    printf("iterNr := %d\n", *iterNr);
    printf("nr of function calls := %d\n", _powell_func_calls);
  }
 
  if(isnan(*fret) || !isfinite(*fret)) {
    if(verbose>10) printf("powell() fails and returns the initial point.\n");
    if(verbose>11) {
      if(isnan(*fret)) printf("  fret := NaN\n");
      if(isfinite(*fret)) printf("  fret := overflow/underflow\n");
    }
    // if failed, then return initial guess
    for(j=0; j<parNr; j++) p[j]=origp[j];
    // and call function again so that any data saved there is correct
    *fret=(*_powellFunc)(parNr, p, _powellFuncData);
    return(3);
  }
  // and call function again so that any data saved there is correct
  *fret=(*_powellFunc)(parNr, p, _powellFuncData);
  if((*iterNr)>=iterMax) return(1);
  if(verbose>0) printf("out of powell() in good order.\n");
  return(0);
}

Referenced by bootstrap(), and tgo().

Variable Documentation

POWELL_LINMIN_MAXIT

int POWELL_LINMIN_MAXIT =100

Max iterations for linear minimization inside powell()

Definition at line 11 of file powell.c.

Referenced by powell(), and tgo().