nlopt.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 "tpcrand.h"
/*****************************************************************************/
#include "tpcnlopt.h"
/*****************************************************************************/
 
/*****************************************************************************/
void nloptInit(
  NLOPT *nlo
) {
  if(nlo==NULL) return;
  nlo->totalNr=0;
  nlo->xfull=NULL;
  nlo->xlower=NULL;
  nlo->xupper=NULL;
  nlo->xdelta=NULL;
  nlo->xtol=NULL;
  nlo->_fun=NULL;
  nlo->fundata=NULL;
  nlo->maxFunCalls=0;
  nlo->funCalls=0;
  nlo->funval=nan("");
  nlo->usePList=0;
  nlo->plist=NULL;
}
/*****************************************************************************/
 
/*****************************************************************************/
void nloptFree(
  NLOPT *nlo
) {
  if(nlo==NULL) return;
  free(nlo->xfull);
  free(nlo->xlower);
  free(nlo->xupper);
  free(nlo->xdelta);
  free(nlo->xtol);
  free(nlo->plist);
  // then set everything to zero or NULL again
  nloptInit(nlo);
}
/*****************************************************************************/
 
/*****************************************************************************/
int nloptAllocate(
  NLOPT *nlo,
  unsigned int parNr
) {
  if(nlo==NULL) return TPCERROR_FAIL;
  /* Delete any previous contents */
  nloptFree(nlo);
  /* If no memory is requested, then just return */
  if(parNr<1) return TPCERROR_OK;
 
  /* Allocate memory for arrays */
  nlo->xfull=calloc(parNr, sizeof(double));
  if(nlo->xfull==NULL) {nloptFree(nlo); return TPCERROR_OUT_OF_MEMORY;}
  nlo->xlower=calloc(parNr, sizeof(double));
  if(nlo->xlower==NULL) {nloptFree(nlo); return TPCERROR_OUT_OF_MEMORY;}
  nlo->xupper=calloc(parNr, sizeof(double));
  if(nlo->xupper==NULL) {nloptFree(nlo); return TPCERROR_OUT_OF_MEMORY;}
  nlo->xdelta=calloc(parNr, sizeof(double));
  if(nlo->xdelta==NULL) {nloptFree(nlo); return TPCERROR_OUT_OF_MEMORY;}
  nlo->xtol=calloc(parNr, sizeof(double));
  if(nlo->xtol==NULL) {nloptFree(nlo); return TPCERROR_OUT_OF_MEMORY;}
  /* Set array length */
  nlo->totalNr=parNr;
  return TPCERROR_OK;
}
/*****************************************************************************/
 
/*****************************************************************************/
int nloptDuplicate(
  NLOPT *nlo1,
  NLOPT *nlo2
) {
  if(nlo1==NULL || nlo2==NULL) return(TPCERROR_FAIL);
  nloptFree(nlo2); if(nlo1->totalNr<1) return(TPCERROR_OK);
 
  int ret=nloptAllocate(nlo2, nlo1->totalNr);
  if(ret!=TPCERROR_OK) return(ret);
 
  nlo2->totalNr=nlo1->totalNr;
  for(unsigned int i=0; i<nlo1->totalNr; i++) nlo2->xfull[i]=nlo1->xfull[i];
  for(unsigned int i=0; i<nlo1->totalNr; i++) nlo2->xlower[i]=nlo1->xlower[i];
  for(unsigned int i=0; i<nlo1->totalNr; i++) nlo2->xupper[i]=nlo1->xupper[i];
  for(unsigned int i=0; i<nlo1->totalNr; i++) nlo2->xdelta[i]=nlo1->xdelta[i];
  for(unsigned int i=0; i<nlo1->totalNr; i++) nlo2->xtol[i]=nlo1->xtol[i];
  nlo2->_fun=nlo1->_fun;
  nlo2->fundata=nlo1->fundata;
  nlo2->maxFunCalls=nlo1->maxFunCalls;
  nlo2->funCalls=nlo1->funCalls;
  nlo2->funval=nlo1->funval;
  nlo2->usePList=nlo1->usePList;
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/
int nloptAddP(
  NLOPT *nlo,
  double *p,
  double funval
) {
  if(nlo==NULL || nlo->totalNr<1 || nlo->funCalls<1) return TPCERROR_FAIL;
  if(nlo->usePList==0) return TPCERROR_OK;
  if(nlo->plist==NULL || nlo->funCalls==1) {
    /* Allocate memory, if this is the first list item */
    nlo->plist=(double*)malloc((nlo->totalNr+1)*sizeof(double));
  } else {
    /* Reallocate memory, if this is not the first list item */
    nlo->plist=(double*)realloc(nlo->plist, (nlo->funCalls)*(nlo->totalNr+1)*sizeof(double));
  }
  if(nlo->plist==NULL) {nloptFree(nlo); return TPCERROR_OUT_OF_MEMORY;}
  /* Copy the values */
  for(unsigned int i=0; i<nlo->totalNr; i++)
    nlo->plist[(nlo->funCalls-1)*(nlo->totalNr+1)+i]=p[i];
  nlo->plist[(nlo->funCalls-1)*(nlo->totalNr+1)+nlo->totalNr]=funval;
  return TPCERROR_OK;
}
/*****************************************************************************/
 
/*****************************************************************************/
int nloptSortP(
  NLOPT *nlo
) {
  if(nlo==NULL) return(TPCERROR_FAIL);
  if(nlo->usePList==0 || nlo->plist==NULL || nlo->totalNr<1) return(TPCERROR_NO_DATA);
  if(nlo->funCalls<2) return(TPCERROR_OK); // nothing to sort
 
  unsigned int dim=nlo->totalNr;
  unsigned int nr=nlo->funCalls;
  for(unsigned int i=0; i<nr-1; i++)
    for(unsigned int j=i+1; j<nr; j++) {
      if( (!isfinite(nlo->plist[i*(dim+1)+dim]) && isfinite(nlo->plist[j*(dim+1)+dim]))
         || nlo->plist[i*(dim+1)+dim] > nlo->plist[j*(dim+1)+dim])
      {
        for(unsigned int k=0; k<=dim; k++) {
          double d=nlo->plist[i*(dim+1)+k];
          nlo->plist[i*(dim+1)+k]=nlo->plist[j*(dim+1)+k]; nlo->plist[j*(dim+1)+k]=d;
        }
      }
    }
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/
int nloptMeanP(
  NLOPT *nlo,
  unsigned int nr,
  double *meanp,
  double *sdp
) {
  if(nlo==NULL || meanp==NULL) return(TPCERROR_FAIL);
  if(nlo->usePList==0 || nlo->plist==NULL || nlo->totalNr<1) return(TPCERROR_NO_DATA);
  if(nlo->funCalls<1) return(TPCERROR_OK);
  if(nr<1 || nr>nlo->funCalls) nr=nlo->funCalls;
 
  unsigned int dim=nlo->totalNr;
  for(unsigned int j=0; j<dim; j++) {
    unsigned int n=0;
    meanp[j]=0.0;
    for(unsigned int i=0; i<nr; i++) {
      if(isfinite(nlo->plist[i*(dim+1)+j])) meanp[j]+=nlo->plist[i*(dim+1)+j];
      n++;
    }
    if(n<1) return(TPCERROR_NO_DATA);
    meanp[j]/=(double)n;
  }
 
  /* Calculate SDs if required */
  if(sdp==NULL) return(TPCERROR_OK);
  for(unsigned int j=0; j<dim; j++) {
    sdp[j]=0.0;
    unsigned int n=0;
    double sqrsum=0.0, sumsqr=0.0;
    for(unsigned int i=0; i<nr; i++) {
      if(isfinite(nlo->plist[i*(dim+1)+j])) {
        sqrsum+=nlo->plist[i*(dim+1)+j];
        sumsqr+=nlo->plist[i*(dim+1)+j]*nlo->plist[i*(dim+1)+j];
        n++;
      }
    }
    if(n<2) continue; // SD=0 if n==1
    sqrsum*=sqrsum;
    double ff=sumsqr - sqrsum/(double)n;
    if(!(ff>0.0)) sdp[j]=0.0;
    else sdp[j]=sqrt( ff / (double)(n-1) );
  }
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/
void nloptPrintP(
  NLOPT *nlo,
  unsigned int nr,
  FILE *fp
) {
  if(nlo==NULL || fp==NULL) return;
  if(nlo->usePList==0 || nlo->plist==NULL || nlo->totalNr<1) return;
  fprintf(fp, "\nSampled points:\n");
  if(nlo->funCalls<1) return;
  if(nr<1 || nr>nlo->funCalls) nr=nlo->funCalls;
 
  unsigned int dim=nlo->totalNr;
  for(unsigned int si=0; si<nr; si++) {
    fprintf(fp, "%d\t", 1+si);
    for(unsigned int i=0; i<dim; i++) fprintf(fp, "%e ", nlo->plist[si*(dim+1)+i]);
    fprintf(fp, "=> %e\n", nlo->plist[si*(dim+1)+dim]);
  }
  fflush(fp);
}
/*****************************************************************************/
 
/*****************************************************************************/
void nloptWrite(
  NLOPT *d,
  FILE *fp
) {
  if(d==NULL || fp==NULL) return;
  if(d->totalNr==0) {fprintf(fp, "NLOPT is empty\n"); return;}
  fprintf(fp, "param      xfull     xlower     xupper     xdelta       xtol\n");
  for(unsigned int i=0; i<d->totalNr; i++) {
    fprintf(fp, "%5d", 1+i);
    fprintf(fp, " %10.5f", d->xfull[i]);
    fprintf(fp, " %10.5f", d->xlower[i]);
    fprintf(fp, " %10.5f", d->xupper[i]);
    fprintf(fp, " %10.5f", d->xdelta[i]);
    fprintf(fp, " %10.5f\n", d->xtol[i]);
  }
  fprintf(fp, "maxFunCalls := %d\n", d->maxFunCalls);
  fprintf(fp, "funCalls := %d\n", d->funCalls);
  fprintf(fp, "funval := %g\n", d->funval);
  fprintf(fp, "usePList := %d\n", d->usePList);
  fflush(fp);
  return;
}
/*****************************************************************************/
 
/*****************************************************************************/
unsigned int nloptLimitFixedNr(
  NLOPT *d
) {
  if(d==NULL || d->totalNr<1) return(0);
  if(d->xlower==NULL || d->xupper==NULL) return(0);
  unsigned int n=0;
  double macheps=doubleMachEps();
  for(unsigned int i=0; i<d->totalNr; i++) {
    double r=d->xupper[i]-d->xlower[i];
    if(!isnan(r) && r<macheps) n++; 
  }
  return(n);
}
/*****************************************************************************/
 
/*****************************************************************************/
unsigned int nloptFixedNr(
  NLOPT *d
) {
  if(d==NULL || d->totalNr<1) return(0);
  unsigned int i, ln=0, dn=0;
  double macheps=doubleMachEps();
  /* based on limits */
  if(d->xlower!=NULL && d->xupper!=NULL) {
    double r;
    for(i=0; i<d->totalNr; i++) {
      r=d->xupper[i]-d->xlower[i];
      if(!isnan(r) && r<macheps) ln++; 
    }
  }
  /* based on deltas */
  if(d->xdelta!=NULL) {
    for(i=0; i<d->totalNr; i++) {
//printf(" %g %g\n", d->xdelta[i], macheps);
      if(!isnan(d->xdelta[i]) && fabs(d->xdelta[i])<macheps) dn++; 
    }
  }
//printf(" %u %u\n", ln, dn);
  /* the one that is nonzero is probably used */
  if(dn==0) return(ln);
  if(ln==0) return(dn);
  /* Both seem to be used; thus check both together */
  unsigned int n=0;
  double r;
  for(i=0; i<d->totalNr; i++) {
    if(!isnan(d->xdelta[i]) && fabs(d->xdelta[i])<macheps) {n++; continue;} 
    r=d->xupper[i]-d->xlower[i];
    if(!isnan(r) && r<macheps) {n++; continue;} 
  }
  return(n);
}
/*****************************************************************************/
 
/*****************************************************************************/
void nloptRemoveEmpties(
  NLOPT *d
) {
  if(d==NULL || d->totalNr<1) return;
  unsigned int i=0, j;
  while(i<d->totalNr) {
    if(!isnan(d->xfull[i]) && !isnan(d->xdelta[i])) {i++; continue;}
    for(j=i+1; j<d->totalNr; j++) {
      d->xfull[j-1]=d->xfull[j];
      d->xlower[j-1]=d->xlower[j];
      d->xupper[j-1]=d->xupper[j];
      d->xdelta[j-1]=d->xdelta[j];
      d->xtol[j-1]=d->xtol[j];
    }
    d->totalNr--;
  }
  return;
}
/*****************************************************************************/
 
/*****************************************************************************/
void nloptdataInit(
  NLOPT_DATA *d
) {
  if(d==NULL) return;
  d->n=0;
  d->x=d->y=d->w=d->sy=NULL;
  d->verbose=0;
  d->fp=stdout;
}
/*****************************************************************************/
 
/*****************************************************************************/