fit_pbr.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"
/*****************************************************************************/
 
/*****************************************************************************/
/* Local functions */
double func_pbr_src(int parNr, double *p, void*);
double func_pbr_fm2(int parNr, double *p, void*);
double func_pbr_rf(int parNr, double *p, void*);
double func_pbr_hill(int parNr, double *p, void*);
/*****************************************************************************/
typedef struct FITDATA {
  unsigned int n;
  double       *x;
  double       *ymeas;
  double       *ysim;
  double       *w;
  optimality_criterion optcrit;
  int           verbose;
} FITDATA;
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Non-linear fitting of a function to sampled plasma-to-blood ratio (PBR) data.",
  " ",
  "Plasma-to-blood ratio is calculated from RBC-to-plasma ratio (RPC) as",
  "  PBR(t)=1/((1-HCT)+HCT*RPC(t))",
  "and RBC-to-plasma ratio is modelled with function 1 (SRC)",
  "  RPC(t)=p1*t*exp(-p2*t) + p3*p1*(1-exp(-p2*t)*(1+p2*t))/(p2*p2)",
  "function 2 (FM2)",
  "  RPC(t)=(p1*t-p3-p5)*exp(-p2*t) + p3*exp(-p4*t) + p5*exp(-p6*t)",
  "function 3 (RF)",
  "  RPC(t)=(p1*t + p3*t^2 + p5*t^3) / (1 + p2*t + p4*t^2 + p6*t^3)",
  "or function 4 (HILL)",
  "  RPC(t)=(p1*t^p2)/(t^p2 + p3)",
  " ",
  "Usage: @P [Options] datafile [parfile]",
  " ",
  "Options:",
  " -HCT=<Haematocrit>",
  "     Constrain haematocrit to specified value; fitted by default.",
  " -model=<SRC|FM2|RF|HILL>",
  "     Use function 1 (SCR, default), function 2 (FM2), function 3 (RF), or",
  "     function 4 (HILL).",
  " -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.",
  " -k=<value>",
  "     Parameter k is constrained to given value.",
  " -lim[=<filename>]",
  "     Specify the constraints for function parameters;",
  "     This file with default values can be created by giving this option",
  "     as the only command-line argument to this program.",
  "     Without file name the default values are printed on screen.",
  " -svg=<Filename>",
  "     Fitted and measured TACs are plotted in specified SVG file.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Datafile must contain at least 2 columns, sample times and ratios,",
  "possibly also weights as last column.",
  " ",
  "Example:",
  "     @P -svg=ia765pbrat.svg ia765pb.rat ia765pbrat.fit",
  " ",
  "See also: fit_bpr, tacinv, tacblend, b2rbc, bpr2cpr, fit_sigm, fit_hiad",
  " ",
  "Keywords: curve fitting, input, blood, plasma, modelling, simulation",
  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    pbrfile[FILENAME_MAX], parfile[FILENAME_MAX];
  char    limfile[FILENAME_MAX], svgfile[FILENAME_MAX];
  int     weights=0; // 0=default, 1=no weighting, 2=frequency
  double  fixed_hct=nan("");
  double  fixed_k=nan("");
  optimality_criterion optcrit=OPTCRIT_OLS;
  unsigned int model=modelCodeIndex("p2bsrc");
 
 
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  pbrfile[0]=parfile[0]=svgfile[0]=limfile[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, "LIM=", 4)==0 && strlen(cptr)>4) {
      strlcpy(limfile, cptr+4, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "SVG=", 4)==0 && strlen(cptr)>4) {
      strlcpy(svgfile, cptr+4, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "K=", 2)==0 && strlen(cptr)>2) {
      if(!atofCheck(cptr+2, &fixed_k) && fixed_k>=0.0) continue;
    } else if(strncasecmp(cptr, "MODEL=", 6)==0 && strlen(cptr)>6) {
      cptr+=6;
      if(strcasecmp(cptr, "SRC")==0) {model=modelCodeIndex("p2bsrc"); continue;}
      if(strcasecmp(cptr, "FM2")==0) {model=modelCodeIndex("p2bfm2"); continue;}
      if(strcasecmp(cptr, "RF")==0) {model=modelCodeIndex("p2brf"); continue;}
      if(strcasecmp(cptr, "HILL")==0) {model=modelCodeIndex("p2bhill"); continue;}
      if(strcasecmp(cptr, "1")==0) {model=modelCodeIndex("p2bsrc"); continue;}
      if(strcasecmp(cptr, "2")==0) {model=modelCodeIndex("p2bfm2"); continue;}
      if(strcasecmp(cptr, "3")==0) {model=modelCodeIndex("p2brf"); continue;}
      if(strcasecmp(cptr, "4")==0) {model=modelCodeIndex("p2bhill"); continue;}
    } else if(strncasecmp(cptr, "HCT=", 4)==0 && strlen(cptr)>4) {
      if(!atofCheck(cptr+4, &fixed_hct) && fixed_hct>0.0 && fixed_hct<1.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(pbrfile, 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);}
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("pbrfile := %s\n", pbrfile);
    if(parfile[0]) printf("parfile := %s\n", parfile);
    if(limfile[0]) printf("limfile := %s\n", limfile);
    if(svgfile[0]) printf("svgfile := %s\n", svgfile);
    printf("weights := %d\n", weights);
    if(!isnan(fixed_k)) printf("fixed_k := %g\n", fixed_k);
    if(!isnan(fixed_hct)) printf("fixed_hct := %g\n", fixed_hct);
    printf("optcrit := %s\n", optcritCode(optcrit));
    printf("model := %s\n", modelDesc(model));
    fflush(stdout);
  }
 
  /* Set parameter names and initial constraints */
  PAR plim; parInit(&plim); if(parAllocate(&plim, modelParNr(model), 1)!=TPCERROR_OK) {
    fprintf(stderr, "Error: cannot initiate parameter list.\n"); return(1);}
  plim.parNr=modelParNr(model);
  strcpy(plim.n[0].name, "HCT"); plim.n[0].lim1=0.36;     plim.n[0].lim2=0.50;
  if(model==modelCodeIndex("p2bsrc")) {
    strcpy(plim.n[1].name, "A");   plim.n[1].lim1=1.0E-06;  plim.n[1].lim2=1.0E+02;
    strcpy(plim.n[2].name, "L");   plim.n[2].lim1=1.0E-06;  plim.n[2].lim2=5.0;
    strcpy(plim.n[3].name, "k");   plim.n[3].lim1=0.0;      plim.n[3].lim2=2.0;
  } else if(model==modelCodeIndex("p2brf")) {
    strcpy(plim.n[1].name, "p1");   plim.n[1].lim1=0.0;  plim.n[1].lim2=2.0;
    strcpy(plim.n[2].name, "p2");   plim.n[2].lim1=0.0;  plim.n[2].lim2=0.5;
    strcpy(plim.n[3].name, "p3");   plim.n[3].lim1=0.0;  plim.n[3].lim2=0.5;
    strcpy(plim.n[4].name, "p4");   plim.n[4].lim1=0.0;  plim.n[4].lim2=1.0;
    strcpy(plim.n[5].name, "p5");   plim.n[5].lim1=0.0;  plim.n[5].lim2=0.1;
    strcpy(plim.n[6].name, "p6");   plim.n[6].lim1=0.0;  plim.n[6].lim2=0.0;
  } else if(model==modelCodeIndex("p2bfm2")) {
    strcpy(plim.n[1].name, "A1");   plim.n[1].lim1=0.0;   plim.n[1].lim2=10.0;
    strcpy(plim.n[2].name, "L1");   plim.n[2].lim1=0.0;   plim.n[2].lim2=5.0;
    strcpy(plim.n[3].name, "A2");   plim.n[3].lim1=0.0;   plim.n[3].lim2=10.0;
    strcpy(plim.n[4].name, "L2");   plim.n[4].lim1=0.0;   plim.n[4].lim2=0.5;
    strcpy(plim.n[5].name, "A3");   plim.n[5].lim1=0.0;   plim.n[5].lim2=5.0;
    strcpy(plim.n[6].name, "L3");   plim.n[6].lim1=-0.1;  plim.n[6].lim2=0.01;
  } else { // if(model==modelCodeIndex("p2bhill"))
    strcpy(plim.n[1].name, "p1");   plim.n[1].lim1=0.0;   plim.n[1].lim2=20.0;
    strcpy(plim.n[2].name, "p2");   plim.n[2].lim1=0.0;   plim.n[2].lim2=5.0;
    strcpy(plim.n[3].name, "p3");   plim.n[3].lim1=0.001; plim.n[3].lim2=5000.0;
  }
  /* If only file name for parameter constraints was given, then write one and exit */
  if(limfile[0] && !pbrfile[0]) {
    /* Check that initial value file does not exist */
    if(strcasecmp(limfile, "stdout")!=0 && fileExist(limfile)) {
      fprintf(stderr, "Error: parameter constraint file %s exists.\n", limfile);
      parFree(&plim); return(1);
    }
    /* Create parameter file */
    plim.parNr=modelParNr(model);
    if(verbose>1 && strcasecmp(limfile, "stdout")!=0) 
      printf("writing parameter constraints file\n");
    if(parWriteLimits(&plim, limfile, verbose-2)!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot write constraints file.\n");
      parFree(&plim); return(1);
    }
    if(strcasecmp(limfile, "stdout")!=0)
      fprintf(stdout, "Parameter constraints written in %s\n", limfile);
    parFree(&plim); return(0);
  }
 
  /* Did we get all the information that we need? */
  if(!pbrfile[0]) { // note that parameter file is optional
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    parFree(&plim); return(1);
  }
  if(optcrit!=OPTCRIT_OLS && optcrit!=OPTCRIT_LAD) {
    fprintf(stderr, "Error: optimality criterion not supported.\n");
    parFree(&plim); return(1);
  }
  if(strcasecmp(limfile, "stdout")==0) limfile[0]=(char)0;
 
 
  /*
   *  Read parameter constraints if file for that was given
   */
  if(limfile[0]) {
    if(parReadLimits(&plim, limfile, verbose-2)!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot read constraints file.\n");
      parFree(&plim); return(1);
    }
    if(verbose>1) parListLimits(&plim, stdout);
  }
  /* Set optional user-defined constraints */
  if(!isnan(fixed_hct)) plim.n[0].lim1=plim.n[0].lim2=fixed_hct;
  if(!isnan(fixed_k) && model==modelCodeIndex("p2bsrc")) plim.n[3].lim1=plim.n[3].lim2=fixed_k;
  if(verbose>1) parListLimits(&plim, stdout);
 
 
 
  /*
   *  Read TAC data
   */
  if(verbose>1) printf("reading %s\n", pbrfile);
  TAC tac; tacInit(&tac);
  if(tacRead(&tac, pbrfile, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); parFree(&plim); return(2);
  }
  if(verbose>1) {
    printf("fileformat := %s\n", tacFormattxt(tac.format));
    printf("tacNr := %d\n", tac.tacNr);
    printf("sampleNr := %d\n", tac.sampleNr);
    printf("xunit := %s\n", unitName(tac.tunit));
    printf("yunit := %s\n", unitName(tac.cunit));
    printf("isframe := %d\n", tac.isframe);
  }
  if(tac.tacNr<1) {
    fprintf(stderr, "Error: file contains no data.\n");
    tacFree(&tac); parFree(&plim); return(2);
  } else if(tac.tacNr>1) {
    if(verbose>0) fprintf(stderr, "Warning: file contains more than one curve.\n");
    tac.tacNr=1;
  }
  if(tac.sampleNr<3) {
    fprintf(stderr, "Error: too few samples.\n");
    tacFree(&tac); parFree(&plim); return(2);
  }
  /* Check NaNs */
  if(tacNaNs(&tac)>0) {
    fprintf(stderr, "Error: data contains missing values.\n");
    tacFree(&tac); parFree(&plim); return(2);
  }
  /* Sort data by sample time */
  tacSortByTime(&tac, &status);
  /* Get x range */
  double xmin, xmax;
  if(tacXRange(&tac, &xmin, &xmax)!=0) {
    fprintf(stderr, "Error: invalid data sample times.\n");
    tacFree(&tac); parFree(&plim); return(2);
  }
  if(verbose>1) {
    printf("xmin := %g\n", xmin);
    printf("xmax := %g\n", xmax);
  }
 
  /* Check and set weights */
  if(weights==0) {
    if(!tacIsWeighted(&tac)) {
      tac.weighting=WEIGHTING_OFF;
      for(int i=0; i<tac.sampleNr; i++) tac.w[i]=1.0;
    } 
  } else if(weights==1) {
    tac.weighting=WEIGHTING_OFF;
    for(int i=0; i<tac.sampleNr; i++) tac.w[i]=1.0;
  } else if(weights==2) {
    if(tacWByFreq(&tac, ISOTOPE_UNKNOWN, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); parFree(&plim); return(2);
    }
  }
 
 
  /*
   *  Copy function parameters into PAR structure for printing and saving
   */
  if(verbose>1) printf("preparing space for parameters\n");
  PAR par; parInit(&par);
  if(parAllocateWithTAC(&par, &tac, modelParNr(model), &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); parFree(&plim); return(3);
  }
  iftFree(&par.h); // remove stupid header info
  /* set time and program name */
  {
    char buf[256];
    tpcProgramName(argv[0], 1, 1, buf, 256);
    iftPut(&par.h, "program", buf, 0, NULL);
    time_t t=time(NULL);
    iftPut(&par.h, "analysis_time", ctime_r_int(&t, buf), 0, NULL);
  }
  par.tacNr=tac.tacNr; par.parNr=modelParNr(model);
  par.format=PAR_FORMAT_TSV_UK;
  for(int i=0; i<par.tacNr; i++) {
    par.r[i].model=model;
    par.r[i].dataNr=tacWSampleNr(&tac);
    par.r[i].start=xmin; 
    par.r[i].end=xmax;
  }
  /* Set parameter names */
  for(int i=0; i<par.parNr; i++) strcpy(par.n[i].name, plim.n[i].name);
  /* set file names */
  iftPut(&par.h, "datafile", pbrfile, 0, NULL);
  /* copy study number */
  {
    int i=iftFindKey(&tac.h, "studynr", 0);
    if(i<0) i=iftFindKey(&tac.h, "study_number", 0);
    if(i>=0) iftPut(&par.h, tac.h.item[i].key, tac.h.item[i].value, 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);
 
  {
    /* Set data pointers for the fit */
    double yfit[tac.sampleNr];
    FITDATA fitdata;
    fitdata.n=tac.sampleNr;
    fitdata.x=tac.x;
    fitdata.ymeas=tac.c[0].y;
    fitdata.ysim=yfit;
    fitdata.w=tac.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 allocate memory.\n");
      tacFree(&tac); parFree(&plim); parFree(&par); return(5);
    }
    if(model==modelCodeIndex("p2bsrc")) nlo._fun=func_pbr_src; 
    else if(model==modelCodeIndex("p2brf")) nlo._fun=func_pbr_rf;
    else if(model==modelCodeIndex("p2bhill")) nlo._fun=func_pbr_hill;
    else if(model==modelCodeIndex("p2bfm2")) nlo._fun=func_pbr_fm2;
    nlo.fundata=&fitdata;
    nlo.totalNr=par.parNr;
    for(int i=0; i<par.parNr; i++) {
      nlo.xlower[i]=plim.n[i].lim1;
      nlo.xupper[i]=plim.n[i].lim2;
    }
    for(int i=0; i<par.parNr; i++) {
      if(!(nlo.xupper[i]>nlo.xlower[i])) nlo.xtol[i]=0.0;
      else nlo.xtol[i]=0.00001*(nlo.xupper[i]-nlo.xlower[i]);
    }
    nlo.maxFunCalls=500000;
    /* Fit */
    if(verbose>4) nloptWrite(&nlo, stdout);
#if(0)
    if(nloptITGO1(&nlo, 1, 0, 10000, 20, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); parFree(&plim); parFree(&par); nloptFree(&nlo); return(6);
    }
#else
    if(nloptIATGO(&nlo, 1, 0, 0.15, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error)); fflush(stderr);
      tacFree(&tac); parFree(&plim); parFree(&par); nloptFree(&nlo); return(6);
    }
#endif
    nlo._fun(nlo.totalNr, nlo.xfull, nlo.fundata);
    if(verbose>2) nloptWrite(&nlo, stdout);
    if(verbose>6) {
      printf("measured and fitted TAC:\n");
      for(int i=0; i<tac.sampleNr; i++)
        printf("\t%g\t%g\t%g\n", tac.x[i], tac.c[0].y[i], yfit[i]);
    }
    par.r[0].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"); fflush(stderr);
      tacFree(&tac); parFree(&plim); 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)); fflush(stderr);
      tacFree(&tac); parFree(&plim); 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");
    /* Prepare place for function TAC */
    double sdist;
    tacGetSampleInterval(&tac, 1.0E-03, &sdist, NULL);
    int snr=(int)simSamples(0.2*sdist, 0.0, par.r[0].end, 2, NULL);
    TAC ftac; tacInit(&ftac);
    int ret=tacAllocate(&ftac, snr, tac.tacNr);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error)); fflush(stderr);
      tacFree(&tac); parFree(&plim); parFree(&par); return(21);
    }
    ftac.tacNr=tac.tacNr; ftac.sampleNr=snr;
    ret=tacCopyHdr(&tac, &ftac);
    for(int ci=0; ci<tac.tacNr && ret==0; ci++)
      ret=tacCopyTacchdr(tac.c+ci, ftac.c+ci);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error)); fflush(stderr);
      tacFree(&tac); parFree(&plim); parFree(&par); tacFree(&ftac);
      return(22);
    }
    snr=(int)simSamples(0.2*sdist, 0.0, par.r[0].end, 2, ftac.x);
    ftac.isframe=0;
    /* Compute function values at sample times */
    ret=0;
    for(int ci=0; ci<tac.tacNr && ret==0; ci++) {
      ret=mfEvalY(modelCode(par.r[ci].model), par.parNr, par.r[ci].p, 
                  ftac.sampleNr, ftac.x, ftac.c[ci].y, verbose-5);
    }
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot calculate function values.\n"); fflush(stderr);
      tacFree(&tac); parFree(&plim); parFree(&par); tacFree(&ftac);
      return(23);
    }
    if(verbose>10) tacWrite(&ftac, stdout, TAC_FORMAT_UNKNOWN, 1, NULL);
    /* Plot */
    ret=tacPlotFitSVG(&tac, &ftac, "", nan(""), nan(""), nan(""), nan(""), svgfile, &status);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error)); fflush(stderr);
      tacFree(&tac); parFree(&plim); parFree(&par); tacFree(&ftac);
      return(24);
    }
    tacFree(&ftac);
    if(verbose>0) printf("Measured and fitted data plotted in %s\n", svgfile);
  }
 
 
 
  tacFree(&tac); parFree(&plim); parFree(&par);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************
 *
 *  Functions to be minimized
 *
 *****************************************************************************/
double func_pbr_src(int parNr, double *p, void *fdata)
{
  FITDATA *d=(FITDATA*)fdata;
 
  if(d->verbose>0) {printf("%s()\n", __func__); fflush(stdout);}
  if(parNr!=4) return(nan(""));
  if(d->verbose>20) printf("p[]: %g %g %g %g\n", p[0], p[1], p[2], p[3]);
 
  /* Calculate the curve values and weighted SS */
  double wss=0.0;
  for(unsigned int i=0; i<d->n; i++) {
    d->ysim[i]=nan("");
    if(isnan(d->ymeas[i]) || isnan(d->x[i]) || d->x[i]<0.0) continue;
    double et=exp(-p[2]*d->x[i]);
    double rc=(1.0-p[3]/p[2])*d->x[i]*et + (p[3]/(p[2]*p[2]))*(1.0-et);
    rc*=p[1];
    d->ysim[i]=1.0/((1.0-p[0])+p[0]*rc);
    double v=d->ysim[i]-d->ymeas[i];
    if(d->optcrit==OPTCRIT_LAD)
      wss+=d->w[i]*fabs(v);
    else
      wss+=d->w[i]*v*v;
  }
 
  return(wss);
}
/*****************************************************************************/
double func_pbr_fm2(int parNr, double *p, void *fdata)
{
  FITDATA *d=(FITDATA*)fdata;
 
  if(d->verbose>0) {printf("%s()\n", __func__); fflush(stdout);}
  if(parNr!=7) return(nan(""));
  if(d->verbose>20) printf("p[]: %g %g %g %g %g %g %g\n", p[0], p[1], p[2], p[3], p[4], p[5], p[6]);
 
  /* Check that D(0)>0 */
  double d0=p[1] - p[3]*p[4] - p[5]*p[6] + p[2]*(p[3]+p[5]);
  if(d0<=0.0) return(nan(""));
 
  /* Calculate the curve values and weighted SS */
  double wss=0.0;
  for(unsigned int i=0; i<d->n; i++) {
    d->ysim[i]=nan("");
    if(!isfinite(d->ymeas[i]) || !isfinite(d->x[i]) || d->x[i]<0.0) continue;
    double rc=(p[1]*d->x[i]-p[3]-p[5])*exp(-p[2]*d->x[i])
              + p[3]*exp(-p[4]*d->x[i]) + p[5]*exp(-p[6]*d->x[i]);
    d->ysim[i]=1.0/((1.0-p[0])+p[0]*rc);
    double v=d->ysim[i]-d->ymeas[i];
    if(d->optcrit==OPTCRIT_LAD)
      wss+=d->w[i]*fabs(v);
    else
      wss+=d->w[i]*v*v;
  }
 
  return(wss);
}
/*****************************************************************************/
double func_pbr_rf(int parNr, double *p, void *fdata)
{
  FITDATA *d=(FITDATA*)fdata;
 
  if(d->verbose>0) {printf("%s()\n", __func__); fflush(stdout);}
  if(parNr!=7) return(nan(""));
  if(d->verbose>20) printf("p[]: %g %g %g %g %g %g %g\n", p[0], p[1], p[2], p[3], p[4], p[5], p[6]);
 
  /* Calculate the curve values and weighted SS */
  double wss=0.0;
  for(unsigned int i=0; i<d->n; i++) {
    d->ysim[i]=nan("");
    if(isnan(d->ymeas[i]) || isnan(d->x[i]) || d->x[i]<0.0) continue;
    double a =       p[1]*d->x[i] + p[3]*d->x[i]*d->x[i] + p[5]*d->x[i]*d->x[i]*d->x[i];
    double b = 1.0 + p[2]*d->x[i] + p[4]*d->x[i]*d->x[i] + p[6]*d->x[i]*d->x[i]*d->x[i];
    d->ysim[i]=1.0/((1.0-p[0])+p[0]*(a/b));
    double v=d->ysim[i]-d->ymeas[i]; if(!isfinite(v)) return(nan(""));
    if(d->optcrit==OPTCRIT_LAD)
      wss+=d->w[i]*fabs(v);
    else
      wss+=d->w[i]*v*v;
  }
 
  return(wss);
}
/*****************************************************************************/
double func_pbr_hill(int parNr, double *p, void *fdata)
{
  FITDATA *d=(FITDATA*)fdata;
 
  if(d->verbose>0) {printf("%s()\n", __func__); fflush(stdout);}
  if(parNr!=4) return(nan(""));
  if(d->verbose>20) printf("p[]: %g %g %g %g\n", p[0], p[1], p[2], p[3]);
 
  /* Calculate the curve values and weighted SS */
  double wss=0.0;
  for(unsigned int i=0; i<d->n; i++) {
    d->ysim[i]=nan("");
    if(isnan(d->ymeas[i]) || isnan(d->x[i]) || d->x[i]<0.0) continue;
    double rc=p[1]*pow(d->x[i], p[2])/(p[3]+pow(d->x[i], p[2]));
    d->ysim[i]=1.0/((1.0-p[0])+p[0]*rc);
    double v=d->ysim[i]-d->ymeas[i]; if(!isfinite(v)) return(nan(""));
    if(d->optcrit==OPTCRIT_LAD)
      wss+=d->w[i]*fabs(v);
    else
      wss+=d->w[i]*v*v;
  }
 
  return(wss);
}
/*****************************************************************************/
 
/*****************************************************************************/