fit_xexp.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
/*****************************************************************************/
#include "tpcextensions.h"
#include "tpcift.h"
#include "tpctac.h"
#include "tpcpar.h"
#include "tpcli.h"
#include "tpccm.h"
#include "tpctacmod.h"
#include "tpclinopt.h"
#include "tpcrand.h"
#include "tpcnlopt.h"
/*****************************************************************************/
 
/*****************************************************************************/
/* Local functions */
double func_xexp(int parNr, double *p, void*);
/*****************************************************************************/
typedef struct FITDATA {
  unsigned int  in;
  double       *ix;
  double       *iy;
  double       *kernel;
  double       *cy;

  unsigned int n;
  double       *x;
  double       *y;
  double       *ysim;
  double       *w;
  int           verbose;
} FITDATA;
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Non-linear fitting of the parameters of response function, convolved with",
  "input TAC, to reproduce given output TAC:",
  "  output(t) = input(t) (x) h(t)",
  "Response function is:",
  "  h(t) = a1*exp(-b1*t) + a2*exp(-b2*t)",
  " ",
  "Usage: @P [Options] inputfile outputfile [parfile]",
  " ",
  "Options:",
  " -a=<value>",
  "     Constrain h(t) integral (AUC) from zero to infinity (a1/b1 + a2/b2)",
  "     to specified value; fitted by default.",
  " -amin=<value>",
  "     Minimum value for AUC; by default 0.1.",
  " -amax=<value>",
  "     Maximum value for AUC; by default 1.0.",
  " -w1",
  "     All weights are set to 1.0 (no weighting); by default, weights in",
  "     output data file are used, if available.",
  " -wf",
  "     Weight by output TAC sampling interval.",
  " -svg=<Filename>",
  "     Fitted and measured TACs are plotted in specified SVG file.",
  " -i=<Interval>",
  "     Sample time interval (sec) in convolution; by default the shortest",
  "     interval in the input data; too long interval as compared to b leads",
  "     to bias.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "If TAC files contain more than one TAC, only the first is used.",
  " ",
  "See also: convexpf, sim_av",
  " ",
  "Keywords: TAC, curve fitting, convolution",
  0};
/*****************************************************************************/
 
/*****************************************************************************/
/* Turn on the globbing of the command line, since it is disabled by default in
   mingw-w64 (_dowildcard=0); in MinGW32 define _CRT_glob instead, if necessary;
   In Unix&Linux wildcard command line processing is enabled by default. */
/*
#undef _CRT_glob
#define _CRT_glob -1
*/
int _dowildcard = -1;
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int     ai, help=0, version=0, verbose=1;
  char    inpfile[FILENAME_MAX], outfile[FILENAME_MAX], parfile[FILENAME_MAX], svgfile[FILENAME_MAX];
  int     weights=0; // 0=default, 1=no weighting, 2=frequency
  double  amin=0.1, amax=1.0;
  double  interval=nan("");
  int     ret;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  inpfile[0]=outfile[0]=parfile[0]=svgfile[0]=(char)0;
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    char *cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(!*cptr) continue;
    if(strcasecmp(cptr, "W1")==0) {
      weights=1; continue;
    } else if(strcasecmp(cptr, "WF")==0) {
      weights=2; continue;
    } else if(strncasecmp(cptr, "SVG=", 4)==0 && strlen(cptr)>4) {
      strlcpy(svgfile, cptr+4, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "AMIN=", 5)==0) {
      amin=atofVerified(cptr+5); if(amin>0.0) continue;
    } else if(strncasecmp(cptr, "AMAX=", 5)==0) {
      amax=atofVerified(cptr+5); if(amax>0.0) continue;
    } else if(strncasecmp(cptr, "A=", 2)==0) {
      amin=amax=atofVerified(cptr+2); if(amin>0.0) continue;
    } else if(strncasecmp(cptr, "I=", 2)==0) {
      interval=atofVerified(cptr+2); if(interval>0.0) continue;
    }
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    return(1);
  } else break;
 
  TPCSTATUS status; statusInit(&status);
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  status.verbose=verbose-3;
  
  /* Print help or version? */
  if(help==2) {tpcHtmlUsage(argv[0], info, ""); return(0);}
  if(help) {tpcPrintUsage(argv[0], info, stdout); return(0);}
  if(version) {tpcPrintBuild(argv[0], stdout); return(0);}
 
  /* Process other arguments, starting from the first non-option */
  if(ai<argc) strlcpy(inpfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(outfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(parfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
  /* Did we get all the information that we need? */
  if(!outfile[0]) { // note that parameter file is optional
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  /* Check parameters */
  if(amin>amax) {
    fprintf(stderr, "Error: invalid limits for a.\n");
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("inpfile := %s\n", inpfile);
    printf("outfile := %s\n", outfile);
    if(parfile[0]) printf("parfile := %s\n", parfile);
    if(svgfile[0]) printf("svgfile := %s\n", svgfile);
    printf("amin := %g\n", amin);
    printf("amax := %g\n", amax);
    printf("weights := %d\n", weights);
    if(!isnan(interval)) printf("interval := %g\n", interval);
    fflush(stdout);
  }
 
 
  /*
   *  Read input and output data
   */
  if(verbose>1) printf("reading TAC data\n");
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  TAC itac, otac; tacInit(&itac); tacInit(&otac);
  int fitSampleNr;
  double fitdur=1.0E+12;
  ret=tacReadModelingData(outfile, inpfile, NULL, NULL, &fitdur, 0,
                          &fitSampleNr, &otac, &itac, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&itac); tacFree(&otac); return(2);
  }
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(itac.format));
    printf("inp.sampleNr := %d\n", itac.sampleNr);
    printf("out.sampleNr := %d\n", otac.sampleNr);
    printf("fitSampleNr := %d\n", fitSampleNr);
    printf("xunit := %s\n", unitName(itac.tunit));
    printf("yunit := %s\n", unitName(itac.cunit));
    printf("fitdur := %g s\n", fitdur);
  }
  if(fitSampleNr<4 || itac.sampleNr<4) {
    fprintf(stderr, "Error: too few samples.\n");
    tacFree(&itac); tacFree(&otac); return(2);
  }
  /* Set second output TAC to place fitted data into */
  otac.tacNr=1;
  ret=tacAllocateMore(&otac, 1);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&itac); tacFree(&otac); return(2);
  }
  /* Check and set weights */
  if(weights==0) {
    if(!tacIsWeighted(&otac)) {
      otac.weighting=WEIGHTING_OFF;
      for(int i=0; i<otac.sampleNr; i++) otac.w[i]=1.0;
    } 
  } else if(weights==1) {
    otac.weighting=WEIGHTING_OFF;
    for(int i=0; i<otac.sampleNr; i++) otac.w[i]=1.0;
  } else if(weights==2) {
    ret=tacWByFreq(&otac, ISOTOPE_UNKNOWN, &status);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&itac); tacFree(&otac); return(2);
    }
  }
 
  /* Sample times should now be in minutes; convert interval to minutes, too. */
  if(interval>0.0) {
    interval/=60.0;
 
    if(interval>0.1*itac.x[itac.sampleNr-1]) {
      fprintf(stderr, "Error: too long interval time.\n");
      tacFree(&itac); tacFree(&otac); return(2);
    }
    if(interval>0.01*itac.x[itac.sampleNr-1]) {
      if(verbose>0) fprintf(stderr, "Warning: interval time may be too long.\n");
    }
  }
 
  /* 
   *  Interpolate data with even sample intervals for convolution
   */
  TAC iitac; tacInit(&iitac);
  if(tacInterpolateToEqualLengthFrames(&itac, interval, interval, &iitac, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: cannot interpolate data to even sample times.\n");
    tacFree(&itac); tacFree(&otac); return(3);
  }
  /* Get the sample interval in interpolated data */
  double freq=iitac.x2[0]-iitac.x1[0]; 
  if(verbose>1) {
    printf("sample_intervals_in_convolution := %g\n", freq);
    printf("interpolated data range: %g - %g\n", iitac.x1[0], iitac.x2[itac.sampleNr-1]);
  }
 
 
  if(verbose>1) fprintf(stdout, "allocate memory for kernel...\n");
  double *kernel=(double*)malloc(2*iitac.sampleNr*sizeof(double));
  if(kernel==NULL) {
    fprintf(stderr, "Error: out of memory.\n");
    tacFree(&itac); tacFree(&otac); tacFree(&iitac); return(4);
  }
  double *cy=kernel+iitac.sampleNr;
 
 
 
 
 
  /*
   *  Fitting
   */
  if(verbose>1) printf("preparing for fitting\n");
 
  /* Set data pointers for the fit */
  FITDATA fitdata;
  fitdata.verbose=0;
  fitdata.kernel=kernel;
  fitdata.in=iitac.sampleNr;
  fitdata.ix=iitac.x;
  fitdata.iy=iitac.c[0].y;
  fitdata.cy=cy;
 
  fitdata.n=otac.sampleNr;
  fitdata.x=otac.x;
  fitdata.y=otac.c[0].y;
  fitdata.ysim=otac.c[1].y;
  fitdata.w=otac.w;
 
 
 
  /* Set parameters for nonlinear optimization */
  /* Because AUC is constrained, we will fit parameters a1, fAUC, a2, and AUC,
     and calculate a1/b1=fAUC*AUC, b1=a1/(fAUC*AUC), and b2=a2/((1-fAUC)*AUC) */
  if(verbose>2) printf("process initial parameter values\n");
  NLOPT nlo; nloptInit(&nlo);
  ret=nloptAllocate(&nlo, 4);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&itac); tacFree(&otac); tacFree(&iitac); return(5);
  }
  /* Set function */
  nlo._fun=func_xexp;
  nlo.fundata=&fitdata;
  nlo.totalNr=4;
  {
    nlo.xlower[0]=1.0E-08; nlo.xupper[0]=amax/freq; nlo.xfull[0]=nlo.xlower[0]+0.02*(nlo.xupper[0]-nlo.xlower[0]);
    nlo.xlower[1]=1.0E-10; nlo.xupper[1]=0.999;   nlo.xfull[1]=0.2;
    nlo.xlower[2]=1.0E-08; nlo.xupper[2]=amax/freq; nlo.xfull[2]=nlo.xlower[2]+0.08*(nlo.xupper[2]-nlo.xlower[2]);
/*
    nlo.xlower[0]=1.0E-08; nlo.xupper[0]=1.0E+00; nlo.xfull[0]=0.6;
    nlo.xlower[1]=1.0E-10; nlo.xupper[1]=0.999;   nlo.xfull[1]=0.2;
    nlo.xlower[2]=1.0E-08; nlo.xupper[2]=1.0E+00; nlo.xfull[2]=0.4;
*/
    nlo.xlower[3]=amin;    nlo.xupper[3]=amax; nlo.xfull[3]=amin+0.8*(amax-amin);
    for(unsigned int i=0; i<nlo.totalNr; i++) {
      nlo.xdelta[i]=0.1*(nlo.xfull[i]-nlo.xlower[i]);
      nlo.xtol[i]=fabs(0.001*nlo.xdelta[i]);
    }
  }
  nlo.maxFunCalls=10000;
 
  /* Get nr of fixed parameters */
  unsigned int fixedNr=nloptFixedNr(&nlo);
  if(verbose>2) printf("fixedNr := %d\n", fixedNr); 
 
#if(0)
  // call function for testing
  fitdata.verbose=10;
  func_xexp(nlo.totalNr, nlo.xfull, &fitdata); exit(1);
#endif
 
  /* Fit */
  drandSeed(1);
  if(verbose>2) nloptWrite(&nlo, stdout);
  if(verbose>0) printf("fitting\n");
  if(nloptSimplexARRS(&nlo, 0, &status)!=TPCERROR_OK) {
  //if(nloptITGO2(&nlo, 1, 0, 0, 0, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&itac); tacFree(&otac); tacFree(&iitac); nloptFree(&nlo); return(6);
  }
  nlo._fun(nlo.totalNr, nlo.xfull, nlo.fundata);
  if(verbose>2) nloptWrite(&nlo, stdout);
  if(verbose>1) {
    printf("measured and fitted TAC:\n");
    for(int i=0; i<otac.sampleNr; i++)
      printf("\t%g\t%g\t%g\n", otac.x[i], otac.c[0].y[i], otac.c[1].y[i]);
  }
  if(verbose>2) printf(" wss := %g\n", nlo.funval);
 
 
  /*
   *  Copy function parameters into PAR structure for printing and saving
   */
  PAR par; parInit(&par);
  ret=parAllocateWithTAC(&par, &otac, nlo.totalNr, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&itac); tacFree(&otac); tacFree(&iitac); nloptFree(&nlo); return(7);
  }
  /* set time and program name */
  {
    char buf[256];
    time_t t=time(NULL);
    iftPut(&par.h, "analysis_time", ctime_r_int(&t, buf), 0, NULL);
    tpcProgramName(argv[0], 1, 1, buf, 256);
    iftPut(&par.h, "program", buf, 0, NULL);
  }
  par.tacNr=1; par.parNr=nlo.totalNr;
  par.format=PAR_FORMAT_TSV_UK;
  par.r[0].model=modelCodeIndex("2exp");
  par.r[0].dataNr=tacWSampleNr(&otac);
  par.r[0].start=otac.x[0]; 
  par.r[0].end=otac.x[otac.sampleNr-1];
  par.r[0].wss=nlo.funval;
  par.r[0].fitNr=nlo.totalNr-fixedNr;
  /* Set parameter names and units */
  for(int i=0; i<par.parNr; i+=2) {
    sprintf(par.n[i].name, "a%d", 1+i/2);
    sprintf(par.n[i+1].name, "b%d", 1+i/2);
    par.n[i].unit=UNIT_UNITLESS;
    par.n[i+1].unit=unitInverse(otac.tunit);
  }
  /* Set parameters */
  par.r[0].p[0]=nlo.xfull[0];
  par.r[0].p[1]=nlo.xfull[0]/(nlo.xfull[1]*nlo.xfull[3]);
  par.r[0].p[2]=nlo.xfull[2];
  par.r[0].p[3]=nlo.xfull[2]/((1.0-nlo.xfull[1])*nlo.xfull[3]);
  // Switch order if necessary to set faster exp first
  if(par.r[0].p[3]>par.r[0].p[1]) {
    double v;
    v=par.r[0].p[0]; par.r[0].p[0]=par.r[0].p[2]; par.r[0].p[2]=v;
    v=par.r[0].p[1]; par.r[0].p[1]=par.r[0].p[3]; par.r[0].p[3]=v;
  }
  /* Set TAC name if missing */
  if(strlen(par.r[0].name)<1) strcpy(par.r[0].name, "1");
  /* set file names */
  iftPut(&par.h, "inputfile", inpfile, 0, NULL);
  iftPut(&par.h, "datafile", outfile, 0, NULL);
 
  /* Print and save */
  if(verbose>0) parWrite(&par, stdout, PAR_FORMAT_TSV_UK, 1, NULL);
  if(parfile[0]) {
    /* Save file */
    if(verbose>1) printf("  saving %s\n", parfile);
    FILE *fp;
    fp=fopen(parfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing.\n");
      tacFree(&itac); tacFree(&otac); tacFree(&iitac); nloptFree(&nlo); parFree(&par); return(11);
    }
    ret=parWrite(&par, fp, PAR_FORMAT_TSV_UK, 1, &status);
    fclose(fp);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&itac); tacFree(&otac); tacFree(&iitac); nloptFree(&nlo); parFree(&par); return(12);
    }
    if(verbose>0) printf("parameters saved in %s\n", parfile);
  }
 
 
 
  /*
   *  Plot measured and fitted data, if requested
   */
  if(svgfile[0]) {
    if(verbose>1) printf("plotting measured and fitted data\n");
    for(unsigned int i=0; i<fitdata.in; i++) iitac.c[0].y[i]=fitdata.cy[i];
    /* Plot */
    ret=tacPlotFitSVG(&otac, &iitac, "", nan(""), nan(""), nan(""), nan(""), svgfile, &status);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&itac); tacFree(&otac); tacFree(&iitac); nloptFree(&nlo); parFree(&par);
      return(24);
    }
 
    if(verbose>0) printf("Measured and fitted data plotted in %s\n", svgfile);
  }
 
 
 
  nloptFree(&nlo); parFree(&par);
  tacFree(&itac); tacFree(&otac); tacFree(&iitac);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************
 *
 *  Functions to be minimized
 *
 *****************************************************************************/
double func_xexp(int parNr, double *p, void *fdata)
{
  FITDATA *d=(FITDATA*)fdata;
 
  if(d->verbose>0) {printf("%s()\n", __func__); fflush(stdout);}
  if(d->verbose>2) printf("p[]: %g %g %g %g\n", p[0], p[1], p[2], p[3]);
 
  /*
   *  Calculate the response function for convolution
   */
  if(d->verbose>1) {fprintf(stdout, "computing the kernel...\n"); fflush(stdout);}
  double a1, b1, a2, b2;
  a1=p[0]; a2=p[2];
  b1=a1/(p[1]*p[3]); b2=a2/((1.0-p[1])*p[3]);
  if(d->verbose>2) printf("a1=%g b1=%g a2=%g b2=%g\n", a1, b1, a2, b2);
  double ksum=0.0;
  double freq=d->ix[1]-d->ix[0];
  if(d->verbose>6) printf("\nData\tKernel:\n");
  for(unsigned int i=0; i<d->in; i++) {
    d->kernel[i] = (a1/b1) * (exp(b1*freq) - 1.0) * exp(-b1*(d->ix[i]+0.5*freq)) +
                   (a2/b2) * (exp(b2*freq) - 1.0) * exp(-b2*(d->ix[i]+0.5*freq));
    ksum+=d->kernel[i];
    if(d->verbose>6) printf("%g\t%g\n", d->ix[i], d->kernel[i]);
  }
  if(d->verbose>2) printf("Sum of response function := %g\n", ksum);
  if(!isnormal(ksum)) {
    if(d->verbose>0) fprintf(stderr, "Error: invalid kernel contents.\n");
    return(nan(""));
  }
 
  /*
   *  Convolution
   */
  if(d->verbose>1) {fprintf(stdout, "convolution...\n"); fflush(stdout);}
  if(convolve1D(d->iy, d->in, d->kernel, d->in, d->cy)!=0) {
    fprintf(stderr, "Error: cannot convolve the data.\n");
    return(nan(""));
  }
  if(d->verbose>4) {
    printf("\nData x\ty\tKernel\tConvolved\n");
    for(unsigned int i=0; i<d->in; i++)
      printf("%g\t%g\t%g\t%g\n", d->ix[i], d->iy[i], d->kernel[i], d->cy[i]);
  }
 
  /* Interpolate to output sample times */
  if(d->verbose>2) {fprintf(stdout, "interpolation...\n"); fflush(stdout);}
  if(liInterpolate(d->ix, d->cy, d->in, d->x, d->ysim, NULL, NULL, d->n, 0, 0, 0)!=0) {
    fprintf(stderr, "Error: cannot interpolate.\n");
    return(nan(""));
  }
  if(d->verbose>3) {
    printf("\nData x\ty\tSimulated\n");
    for(unsigned int i=0; i<d->n; i++)
      printf("%g\t%g\t%g\n", d->x[i], d->y[i], d->ysim[i]);
  }
 
  /* Calculate the weighted SS */
  if(d->verbose>2) {fprintf(stdout, "computing WSS...\n"); fflush(stdout);}
  double wss=0.0;
  for(unsigned i=0; i<d->n; i++) {
    double v=d->ysim[i] - d->y[i];
    wss+=d->w[i]*v*v;
  }
  if(d->verbose>1) {
    for(int i=0; i<parNr; i++) printf(" %g", p[i]);
    printf(" -> %g\n", wss);
  }
 
  return(wss);
}
/*****************************************************************************/
 
/*****************************************************************************/