nlopt1d.c

Go to the documentation of this file.

 
/******************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
/*****************************************************************************/
#include "tpcextensions.h"
/*****************************************************************************/
#include "tpcnlopt.h"
/*****************************************************************************/
 
/*****************************************************************************/
int nlopt1D(
  NLOPT *nlo,
  unsigned int maxIter,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>1) {printf("%s(NLOPT, %d, status)\n", __func__, maxIter); fflush(stdout);}
  if(nlo==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  if(nlo->totalNr<1 || nlo->xfull==NULL || nlo->_fun==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
 
  /* Check the number of fixed parameters, based on constraints or delta */
  if(nlo->totalNr-nloptFixedNr(nlo)!=1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_PARNR);
    return TPCERROR_INVALID_PARNR;
  }
  /* Which parameter is the free one? */
  double macheps=doubleMachEps();
  unsigned int fpi;
  for(fpi=0; fpi<nlo->totalNr; fpi++) {
    double r=nlo->xupper[fpi]-nlo->xlower[fpi];
    if(isfinite(nlo->xdelta[fpi]) && fabs(nlo->xdelta[fpi])>macheps && isfinite(r) && r>macheps)
      break; 
  }
  if(fpi==nlo->totalNr) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_PARNR);
    return TPCERROR_INVALID_PARNR;
  }
  if(verbose>3) printf("  index_of_free_parameter := %u\n", fpi);
 
  /* Check that initial value is inside its limits */
  double pmin=nlo->xlower[fpi];
  double pmax=nlo->xupper[fpi];
  double delta=nlo->xdelta[fpi];
  double tol=nlo->xtol[fpi];
  if(verbose>2)
    printf("  initial_parameter := %g [%g, %g, %g]\n", nlo->xfull[fpi], pmin, pmax, delta);
  if(!(nlo->xfull[fpi]>=pmin) || !(nlo->xfull[fpi]<=pmax)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
  /* Check that limits are wider than 2*delta */
  if(!( (pmax-pmin) > 2.0*delta)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
 
  /* Set the maximum number of iterations, if not given by user */
  if(maxIter<100) maxIter=200+(pmax-pmin)/delta;
  if(verbose>1) printf("  maxIter := %u\n", maxIter);
 
  /* Set 3 bracketing points */
  double p1=0., p2=0., p3=0., f1=0., f2=0., f3=0.;
  double minf, minp;
  /* Set an array function parameters */
  double p[nlo->totalNr];
  for(unsigned int i=0; i<nlo->totalNr; i++) p[i]=nlo->xfull[i];
 
  /* Calculate function value with initial parameter value */
  unsigned int iterNr=1;
  p2=minp=nlo->xfull[fpi]; if(p2<pmin+delta) p2=pmin+delta; else if(p2>pmax-delta) p2=pmax-delta;
  p[fpi]=p2; minf=f2=(*nlo->_fun)(nlo->totalNr, p, nlo->fundata);
  if(!isfinite(f2)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
  nlo->funCalls++;
  if(nlo->usePList) {
    int ret=nloptAddP(nlo, p, f2);
    if(ret!=TPCERROR_OK) {statusSet(status, __func__, __FILE__, __LINE__, ret); return(ret);}
  }
 
  /* Set the two other bracketing points and compute their function values */
  p1=p2-delta; p[fpi]=p1; f1=(*nlo->_fun)(nlo->totalNr, p, nlo->fundata); iterNr++;
  nlo->funCalls++;
  if(nlo->usePList) {
    int ret=nloptAddP(nlo, p, f1);
    if(ret!=TPCERROR_OK) {statusSet(status, __func__, __FILE__, __LINE__, ret); return(ret);}
  }
  if(f1<minf) {minf=f1; minp=p1;}
 
  p3=p2+delta; p[fpi]=p3; f3=(*nlo->_fun)(nlo->totalNr, p, nlo->fundata); iterNr++;
  nlo->funCalls++;
  if(nlo->usePList) {
    int ret=nloptAddP(nlo, p, f3);
    if(ret!=TPCERROR_OK) {statusSet(status, __func__, __FILE__, __LINE__, ret); return(ret);}
  }
  if(f3<minf) {minf=f3; minp=p3;}
 
 
  /* Now we have 3 points on the function. 
     Start looking for a bracketing set such that f1 > f2 < f3 is the case. */
  double jump_size=delta;
  while(!(f1>f2 && f2<f3)) {
    if(verbose>3) {
      printf("  jump_size := %g\n", jump_size);
      printf("  bracketing points: %g %g %g -> %g %g %g\n", p1, p2, p3, f1, f2, f3);
      printf("  min: %g -> %g\n", minp, minf);
    }
    /* check against stopping rules */
    if(iterNr>=maxIter) break;
    if((p3-p1)<tol) break;
    /* if f1 is small then take a step to the left */
    if(f1<f3) { 
      /* check if the minimum is colliding against the bounds. If so then pick a point between 
         p1 and p2 in the hopes that shrinking the interval will be a good thing to do.
         Or if p1 and p2 aren't differentiated then try and get them to obtain different values. */
      if(p1==pmin || (fabs(f1-f2)<macheps && (pmax-pmin)<jump_size )) {
        p3=p2; f3=f2; p2=0.5*(p1+p2); 
        p[fpi]=p2; f2=(*nlo->_fun)(nlo->totalNr, p, nlo->fundata); iterNr++;
        nlo->funCalls++;
        if(nlo->usePList) {
          int ret=nloptAddP(nlo, p, f2);
          if(ret!=TPCERROR_OK) {statusSet(status, __func__, __FILE__, __LINE__, ret); return(ret);}
        }
        if(f2<minf) {minf=f2; minp=p2;}
      } else {
        /* pick a new point to the left of our current bracket */
        p3=p2; f3=f2; p2=p1; f2=f1;
        p1-=jump_size; if(p1<pmin) p1=pmin;
        p[fpi]=p1; f1=(*nlo->_fun)(nlo->totalNr, p, nlo->fundata); iterNr++;
        nlo->funCalls++;
        if(nlo->usePList) {
          int ret=nloptAddP(nlo, p, f1);
          if(ret!=TPCERROR_OK) {statusSet(status, __func__, __FILE__, __LINE__, ret); return(ret);}
        }
        if(f1<minf) {minf=f1; minp=p1;}
        jump_size*=2.0;
      }
    } else { // otherwise f3 is small and we should take a step to the right.
      /* check if the minimum is colliding against the bounds. If so then pick a point between 
         p2 and p3 in the hopes that shrinking the interval will be a good thing to do.
         Or if p2 and p3 aren't differentiated then try and get them to obtain different values. */
      if(p3==pmax || (fabs(f2-f3)<macheps && (pmax-pmin)<jump_size)) {
        p1=p2; f1=f2; p2=0.5*(p3+p2);
        p[fpi]=p2; f2=(*nlo->_fun)(nlo->totalNr, p, nlo->fundata); iterNr++;
        nlo->funCalls++;
        if(nlo->usePList) {
          int ret=nloptAddP(nlo, p, f2);
          if(ret!=TPCERROR_OK) {statusSet(status, __func__, __FILE__, __LINE__, ret); return(ret);}
        }
        if(f2<minf) {minf=f2; minp=p2;}
      } else {
        /* pick a new point to the right of our current bracket */
        p1=p2; f1=f2; p2=p3; f2=f3;
        p3+=jump_size; if(p3>pmax) p3=pmax;
        p[fpi]=p3; f3=(*nlo->_fun)(nlo->totalNr, p, nlo->fundata); iterNr++;
        nlo->funCalls++;
        if(nlo->usePList) {
          int ret=nloptAddP(nlo, p, f3);
          if(ret!=TPCERROR_OK) {statusSet(status, __func__, __FILE__, __LINE__, ret); return(ret);}
        }
        if(f3<minf) {minf=f3; minp=p3;}
        jump_size*=2.0;
      }
    }
  }
  if(verbose>3) {
    printf("  brackets ready.\n\n");
    if(verbose>3) {
      printf("  bracketing points: %g %g %g -> %g %g %g\n", p1, p2, p3, f1, f2, f3);
      printf("  min: %g -> %g\n", minp, minf);
    }
  }
 
  /* Loop until we have done the max allowable number of iterations or the bracketing window is 
     smaller than required tolerance.
     Within this loop we maintain the invariant that: f1 > f2 < f3 and p1 < p2 < p3. */
  const double tau=0.05;
  while(iterNr<maxIter && (p3-p1)>tol) {
 
    //p_min = lagrange_poly_min_extrap(p1,p2,p3, f1,f2,f3);
    double d=f1*(p3*p3-p2*p2) + f2*(p1*p1-p3*p3) + f3*(p2*p2-p1*p1);
    double e=2.0*(f1*(p3-p2) + f2*(p1-p3) + f3*(p2-p1));
    if(fabs(e)<macheps || !isfinite(d/=e)) {
      minp=p2;
    } else {
      if(p1<=d && d<=p3) {minp=d;} else {minp=d; if(p1>minp) minp=p1; if(p3<minp) minp=p3;}
    }
 
    /* make sure min p isn't too close to the three points we already have */
    if(minp<p2) {
      d=(p2-p1)*tau; if(fabs(p1-minp)<d) minp=p1+d; else if(fabs(p2-minp)<d) minp=p2-d;
    } else {
      d=(p3-p2)*tau; if(fabs(p2-minp)<d) minp=p2+d; else if(fabs(p3-minp)<d) minp=p3-d;
    }
 
    /* make sure one side of the bracket isn't super huge compared to the other side.
       If it is then contract it. */
    double bracket_ratio=fabs(p1-p2)/fabs(p2-p3);
    if(!(bracket_ratio<100.0 && bracket_ratio>0.01)) {
      /* Force min p to be on a reasonable side. */
      if(bracket_ratio>1.0 && minp>p2) minp=0.5*(p1+p2); else if(minp<p2) minp=0.5*(p2+p3);
    }
 
    /* Compute function value at min p */
    p[fpi]=minp; minf=(*nlo->_fun)(nlo->totalNr, p, nlo->fundata); iterNr++;
    nlo->funCalls++;
    if(nlo->usePList) {
      int ret=nloptAddP(nlo, p, minf);
      if(ret!=TPCERROR_OK) {statusSet(status, __func__, __FILE__, __LINE__, ret); return(ret);}
    }
 
    /* Remove one of the endpoints of our bracket depending on where the new point falls */
    if(minp<p2) {
      if(f1>minf && minf<f2) {p3=p2; f3=f2; p2=minp; f2=minf;} else {p1=minp; f1=minf;}
    } else {
      if(f2>minf && minf<f3) {p1=p2; f1=f2; p2=minp; f2=minf;} else {p3=minp; f3=minf;}
    }
 
    if(verbose>3) {
      printf("  bracketing points: %g %g %g -> %g %g %g\n", p1, p2, p3, f1, f2, f3);
      printf("  min: %g -> %g\n", minp, minf);
    }
  }
 
  /* Copy the estimated parameter into nlo structure */
  nlo->xfull[fpi]=p2; // Not minp!
  /* Compute function value at the minimum, in case user uses the fitted data in nlo structure */
  nlo->funval=(*nlo->_fun)(nlo->totalNr, nlo->xfull, nlo->fundata);
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return TPCERROR_OK;
}
/*****************************************************************************/
 
/*****************************************************************************/