fit_disp.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 "tpcfileutil.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 "tpcmodels.h"
#include "tpcnlopt.h"
/*****************************************************************************/
 
/*****************************************************************************/
const int maxParNr=2;
/* Local functions */
double func_disp(int parNr, double *p, void*);
/*****************************************************************************/
typedef struct FITDATA {
  unsigned int ni;
  double       *xi;
  double       *yi;

  unsigned int no;
  double       *xo;
  double       *yo;
  double       *w;
  double       *ysim;

  optimality_criterion optcrit;
  int           verbose;
} FITDATA;
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Non-linear fitting of dispersion and delay from system input and output",
  "curves.",
  " ",
  "Usage: @P [Options] inputfile outputfile [parfile]",
  " ",
  "Options:",
  " -min=<OLS|LAD>",
  "     Sum-of-squares (OLS) is minimized by default, but optionally",
  "     sum of absolute deviations (LAD) can be selected.",
  " -w1",
  "     All weights are set to 1.0 (no weighting); by default, weights in",
  "     data file are used, if available.",
  " -wf",
  "     Weight by sampling interval.",
  " -svg=<Filename>",
  "     Fitted and measured TACs are plotted in specified SVG file.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Data files must contain 2 columns, sample times and concentrations,",
  "possibly also weights as last column.",
  " ",
  "Example:",
  "     @P -svg=delayfit.svg input.tac output.tac delayfit.par",
  " ",
  "See also: disp4dft, fit_sinf, fit_xexp, convexpf, fit_wrlv",
  " ",
  "Keywords: dispersion, curve fitting, input, blood",
  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    tacfile1[FILENAME_MAX], tacfile2[FILENAME_MAX], 
          parfile[FILENAME_MAX], svgfile[FILENAME_MAX];
  int     weights=0; // 0=default, 1=no weighting, 2=frequency
  optimality_criterion optcrit=OPTCRIT_OLS;
 
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  tacfile1[0]=tacfile2[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, "MIN=", 4)==0 && strlen(cptr)>4) {
      optcrit=optcritId(cptr+4); if(optcrit!=OPTCRIT_UNKNOWN) continue;
    } else if(strncasecmp(cptr, "SVG=", 4)==0 && strlen(cptr)>4) {
      strlcpy(svgfile, cptr+4, FILENAME_MAX); 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(tacfile1, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(tacfile2, 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(!tacfile2[0]) { // note that parameter file is optional
    fprintf(stderr, "Error: missing file name.\n");
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("tacfile1 := %s\n", tacfile1);
    printf("tacfile2 := %s\n", tacfile2);
    if(parfile[0]) printf("parfile := %s\n", parfile);
    if(svgfile[0]) printf("svgfile := %s\n", svgfile);
    printf("weights := %d\n", weights);
    printf("optcrit := %s\n", optcritCode(optcrit));
    fflush(stdout);
  }
 
 
  /*
   *  Read TAC data
   */
  if(verbose>1) printf("reading TAC data\n");
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  TAC tac1, tac2; tacInit(&tac1); tacInit(&tac2);
  int fitSampleNr;
  double fitdur=1.0E+12;
  if(tacReadModelingData(tacfile2, tacfile1, NULL, NULL, &fitdur, 0,
                          &fitSampleNr, &tac2, &tac1, &status) != TPCERROR_OK)
  {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac1); tacFree(&tac2); return(2);
  }
  /* Convert sample times into seconds */
  if(tacXUnitConvert(&tac1, UNIT_SEC, &status)!=TPCERROR_OK ||
     tacXUnitConvert(&tac2, UNIT_SEC, &status)!=TPCERROR_OK) 
  {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac1); tacFree(&tac2); return(2);
  }
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(tac1.format));
    printf("tac1.sampleNr := %d\n", tac1.sampleNr);
    printf("tac2.sampleNr := %d\n", tac2.sampleNr);
    printf("fitSampleNr := %d\n", fitSampleNr);
    printf("xunit := %s\n", unitName(tac1.tunit));
    printf("yunit := %s\n", unitName(tac1.cunit));
    printf("fitdur := %g s\n", fitdur);
  }
  if(fitSampleNr<8 || tac1.sampleNr<8 || tac2.sampleNr<8) {
    fprintf(stderr, "Error: too few samples.\n");
    tacFree(&tac1); tacFree(&tac2); return(2);
  }
 
 
  /* Check and set weights */
  if(weights==0) {
    if(!tacIsWeighted(&tac2)) {
      tac2.weighting=WEIGHTING_OFF;
      for(int i=0; i<tac2.sampleNr; i++) tac2.w[i]=1.0;
    } 
  } else if(weights==1) {
    tac2.weighting=WEIGHTING_OFF;
    for(int i=0; i<tac2.sampleNr; i++) tac2.w[i]=1.0;
  } else if(weights==2) {
    if(tacWByFreq(&tac2, ISOTOPE_UNKNOWN, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac1); tacFree(&tac2); return(2);
    }
  }
 
 
  /*
   *  Copy function parameters into PAR struct for printing and saving
   */
  if(verbose>1) printf("preparing space for parameters\n");
  PAR par; parInit(&par);
  if(parAllocateWithTAC(&par, &tac2, maxParNr, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac1); tacFree(&tac2); parFree(&par); return(3);
  }
  iftFree(&par.h); // remove stupid header info
  /* 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=2;
  par.format=PAR_FORMAT_TSV_UK;
  for(int i=0; i<par.tacNr; i++) {
    par.r[i].model=modelCodeIndex("dispdelay");
    par.r[i].dataNr=tacWSampleNr(&tac2);
    par.r[i].start=tac2.x[0]; 
    par.r[i].end=tac2.x[tac2.sampleNr-1];
  }
  /* Set parameter names */
  strcpy(par.n[0].name, "tau"); par.n[0].unit=UNIT_SEC;
  strcpy(par.n[1].name, "delay"); par.n[1].unit=UNIT_SEC;
  /* set file names */
  iftPut(&par.h, "inputfile", tacfile1, 0, NULL);
  iftPut(&par.h, "datafile", tacfile2, 0, NULL);
  /* set optimality criterion */
  iftPut(&par.h, "optimality_criterion", optcritCode(optcrit), 0, NULL);
 
 
  /*
   *  Fitting
   */
  if(verbose>1) printf("preparing for fitting\n");
 
  drandSeed(1);
 
  double yfit[tac2.sampleNr];
  FITDATA fitdata;
  {
    /* Set data pointers for the fit */
    fitdata.ni=tac1.sampleNr;
    fitdata.xi=tac1.x;
    fitdata.yi=tac1.c[0].y;
    fitdata.no=tac2.sampleNr;
    fitdata.xo=tac2.x;
    fitdata.yo=tac2.c[0].y;
    fitdata.ysim=yfit;
    fitdata.w=tac2.w;
    fitdata.optcrit=optcrit;
    if(verbose>10) fitdata.verbose=verbose-10; else fitdata.verbose=0;
    /* Set NLLS options */
    NLOPT nlo; nloptInit(&nlo);
    if(nloptAllocate(&nlo, par.parNr)!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot initiate NLLS.\n");
      tacFree(&tac1); tacFree(&tac2); parFree(&par); return(5);
    }
    nlo._fun=func_disp;
    nlo.fundata=&fitdata;
    nlo.totalNr=par.parNr;
    /* Set initial values and limits */
    nlo.xlower[0]=0.0; nlo.xupper[0]=60.0; nlo.xfull[0]=10.0; nlo.xdelta[0]=1.0; nlo.xtol[0]=0.02; 
    nlo.xlower[1]=0.0; nlo.xupper[1]=60.0; nlo.xfull[1]=10.0; nlo.xdelta[1]=1.0; nlo.xtol[1]=0.02; 
 
    if(verbose>2) printf("non-linear optimization\n");
    if(nloptSimplex(&nlo, 0, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac1); tacFree(&tac2); parFree(&par); nloptFree(&nlo); return(6);
    }
    if(verbose>3) {
      nloptWrite(&nlo, stdout);
      (void)func_disp(nlo.totalNr, nlo.xfull, nlo.fundata);
      if(verbose>6) {
        printf("measured and fitted TAC:\n");
        for(unsigned int i=0; i<fitdata.no; i++)
          printf("\t%g\t%g\t%g\t%g\n", fitdata.yi[i], fitdata.yo[i], fitdata.ysim[i], fitdata.w[i]);
      }
      double final_wss=0.0;
      for(unsigned int i=0; i<fitdata.no; i++) {
        double v=fitdata.ysim[i]-fitdata.yo[i];
        final_wss+=fitdata.w[i]*v*v;
      }
      printf(" final_wss := %g\n", final_wss);
    }
 
    /* Simplex again, with new deltas */
    for(unsigned int i=0; i<nlo.totalNr; i++) {
      if(fabs(nlo.xdelta[i])<1.0E-100) continue;
      if(nlo.xfull[i]>1.0E-06) nlo.xdelta[i]=0.0103*nlo.xfull[i];
      else nlo.xdelta[i]*=0.0103;
      nlo.xtol[i]*=0.5;
    }
    if(verbose>2) printf("non-linear optimization, 2nd time\n");
    if(nloptSimplex(&nlo, 0, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac1); tacFree(&tac2); parFree(&par); nloptFree(&nlo); return(6);
    }
    if(verbose>3) {
      nloptWrite(&nlo, stdout);
      (void)func_disp(nlo.totalNr, nlo.xfull, nlo.fundata);
      if(verbose>6) {
        printf("measured and fitted TAC:\n");
        for(unsigned int i=0; i<fitdata.no; i++)
          printf("\t%g\t%g\t%g\t%g\n", fitdata.yi[i], fitdata.yo[i], fitdata.ysim[i], fitdata.w[i]);
      }
      double final_wss=0.0;
      for(unsigned int i=0; i<fitdata.no; i++) {
        double v=fitdata.ysim[i]-fitdata.yo[i];
        final_wss+=fitdata.w[i]*v*v;
      }
      printf(" final_wss := %g\n", final_wss);
    }
 
    /* and Simplex again, with new deltas */
    for(unsigned int i=0; i<nlo.totalNr; i++) {
      nlo.xdelta[i]=0.0093*nlo.xdelta[i];
      nlo.xtol[i]*=0.001;
    }
    if(verbose>2) printf("non-linear optimization, 3rd time\n");
    if(nloptSimplex(&nlo, 0, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac1); tacFree(&tac2); parFree(&par); nloptFree(&nlo); return(6);
    }
 
    /* Call objective function again with the final parameters */
    (void)func_disp(nlo.totalNr, nlo.xfull, nlo.fundata);
    if(verbose>3) nloptWrite(&nlo, stdout);
    if(verbose>6) {
      printf("measured and fitted TAC:\n");
      for(unsigned int i=0; i<fitdata.no; i++)
        printf("\t%g\t%g\n", fitdata.yo[i], fitdata.ysim[i]);
    }
    double final_wss=0.0;
    {
      for(unsigned int i=0; i<fitdata.no; i++) {
        double v=fitdata.ysim[i]-fitdata.yo[i];
        final_wss+=fitdata.w[i]*v*v;
      }
    }
    if(verbose>3) printf(" final_final_final_wss := %g\n", final_wss);
    par.r[0].wss=final_wss; //nlo.funval;
 
    /* Copy parameters */
    for(int i=0; i<par.parNr; i++) par.r[0].p[i]=nlo.xfull[i];
 
    nloptFree(&nlo);
  }
 
  /* Print and save parameters */
  if(verbose>0 || !parfile[0]) parWrite(&par, stdout, PAR_FORMAT_TSV_UK, 1, NULL);
  if(parfile[0]) {
    par.format=parFormatFromExtension(parfile);
    if(verbose>2) printf("parameter file format := %s\n", parFormattxt(par.format));
    if(par.format==PAR_FORMAT_UNKNOWN) par.format=PAR_FORMAT_TSV_UK;
    /* Save file */
    if(verbose>1) printf("  saving %s\n", parfile);
    FILE *fp=fopen(parfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing.\n");
      tacFree(&tac1); tacFree(&tac2); parFree(&par); return(11);
    }
    int ret=parWrite(&par, fp, PAR_FORMAT_UNKNOWN, 1, &status);
    fclose(fp);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac1); tacFree(&tac2); 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");
    TAC fit; tacInit(&fit);
    (void)tacDuplicate(&tac2, &fit);
    for(int i=0; i<tac2.sampleNr; i++) fit.c[0].y[i]=fitdata.ysim[i];
    /* Plot */
    if(tacPlotFitSVG(&tac2, &fit, "", nan(""), nan(""), nan(""), nan(""), svgfile, &status)
       != TPCERROR_OK)
    {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac1); tacFree(&tac2); parFree(&par); tacFree(&fit);
      return(21);
    }
    if(verbose>0) printf("Measured and fitted data plotted in %s\n", svgfile);
    tacFree(&fit);
  }
 
 
  tacFree(&tac1); tacFree(&tac2); parFree(&par);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************
 *
 *  Functions to be minimized
 *
 *****************************************************************************/
double func_disp(int parNr, double *p, void *fdata)
{
  FITDATA *d=(FITDATA*)fdata;
 
  if(d->verbose>0) {printf("%s()\n", __func__); fflush(stdout);}
  if(parNr!=2) return(nan(""));
  if(d->verbose>20) printf("p[]: %g %g\n", p[0], p[1]);
 
  /* Add dispersion to the input TAC */
  double dtac[2*d->ni];
  double *buf=dtac+d->ni;
  for(unsigned int i=0; i<d->ni; i++) dtac[i]=d->yi[i];
  if(simDispersion(d->xi, dtac, d->ni, p[0], 0.0, buf)) return(nan(""));
 
  /* Add delay and interpolate to output sample times */
  for(unsigned int i=0; i<d->ni; i++) buf[i]=d->xi[i]+p[1];
  if(liInterpolate(buf, dtac, d->ni, d->xo, d->ysim, NULL, NULL, d->no, 0, 1, 0)!=0) return(nan(""));
 
  /* Calculate the weighted SS */
  double wss=0.0;
  for(unsigned int i=0; i<d->no; i++) {
    if(isnan(d->ysim[i]) || isnan(d->xo[i]) || d->xo[i]<0.0) continue;
    double v=d->ysim[i]-d->yo[i];
    if(d->optcrit==OPTCRIT_LAD)
      wss+=d->w[i]*fabs(v);
    else
      wss+=d->w[i]*v*v;
  }
 
  return(wss);
}
/*****************************************************************************/
 
/*****************************************************************************/