fit_sigm.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcmodel.h"
#include "libtpccurveio.h"
#include "libtpcsvg.h"
#include "libtpcmodext.h"
/*****************************************************************************/
int FIT2DAT_MAXNR=10000;
/*****************************************************************************/
 
/*****************************************************************************/
double func841(int parNr, double *p, void*);
double func844(int parNr, double *p, void*);
double func2801(int parNr, double *p, void*);
/*****************************************************************************/
double *x, *ymeas, *yfit, *w; /* Local */
int dataNr=0, parNr=3;
double pmin[MAX_PARAMETERS], pmax[MAX_PARAMETERS];
double min1, min2, max1, max2;
double (*func)(int parNr, double *p, void*);
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Non-linear fitting of the Hill function (sigmoid function) to a (xi,yi)",
  "curve data. Optionally, a background term B is added to the function:",
  " ",
  "       f(x) = (A * x^n ) / ( x^n + K ) + B , where A>=0, K>0, n>0, B>=0 ",
  " ",
  "Function for estimating EC50 from sigmoidal dose-response curve is obtained",
  "by rearranging and marking n=HillSlope, A=Top-Bottom, B=Bottom, and",
  "K=EC50^n :",
  " ",
  "       f(x) = Bottom + (Top-Bottom)/(1 + (EC50/x)^HillSlope)",
  " ",
  "With HillSlope constrained to 1 (see options below) the function reduces to",
  "unconstrained one site binding hyperbola, f(x)=Top*x/(EC50+x).",
  " ",
  "Usage: @P [Options] tacfile [fitfile]",
  " ",
  "Options:",
  " -EC50",
  "     EC50 is estimated by fitting sigmoidal dose-response curve function",
  "     when drug concentrations are on a linear scale (not on a log scale);",
  "     By default, parameter Top is fitted and Bottom is constrained to 0,",
  "     which can be changed with options (see below).",
  " -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.",
  " -n=<value>",
  "     Parameter n or HillSlope is constrained to the given value; set -n=1",
  "     to fit one site binding hyperbola.",
  " -B  Parameter B or Bottom is fitted; by default B=0",
  " -B=<value>",
  "     Parameter B or Bottom is constrained to the given value.",
  " -A=<value>",
  "     Parameter Top is constrained to the given value; currently this option",
  "     does not affect parameter A in the default function.",
  " -MRL",
  "     Error is returned if MRL check is not passed.",
  " -res=<Filename>",
  "     Fitted parameters are also written in result file format (3).",
  " -svg=<Filename>",
  "     Fitted and measured TACs are plotted in specified SVG file.",
  " -fit=<Filename>",
  "     Fitted TACs are written in DFT format.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "TAC file must contain at least 2 columns, sample times (or concentrations)",
  "and measurement values (or receptor occupancies).",
  "Weights can be specified as usual if data is in DFT format (1).",
  " ",
  "Program writes the fit start and end times, nr of points, WSS,",
  "and parameters of the fitted function to the fit file (2).",
  " ",
  "References:",
  "1. https://www.turkupetcentre.net/analysis/doc/format_dft.html",
  "2. https://www.turkupetcentre.net/analysis/doc/format_fit.html",
  "3. https://www.turkupetcentre.net/analysis/doc/format_res.html",
  " ",
  "See also: fit_hiad, fit_ppf, fit_exp, fit_ratf, fit_bpr, fit2dat, tacinv",
  " ",
  "Keywords: curve fitting, Hill function, sigmoid, EC50, dose-response curve",
  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;
  int      fi, pi, ri, type=841, m, n, ret;
  int      is_n_fixed=0, is_b_fixed=1, is_a_fixed=0;
  int      weights=0; // 0=default, 1=no weighting, 2=frequency
  double   fixed_n=1.0, fixed_b=0.0, fixed_a=100.0;
  int      MRL_check=0; // 0=off, 1=on
  char     datfile[FILENAME_MAX], fitfile[FILENAME_MAX], outfile[FILENAME_MAX],
           svgfile[FILENAME_MAX], resfile[FILENAME_MAX], *cptr, temp[512];
  char     progname[256];
  DFT      dft;
  FIT      fit;
  RES      res;
  double   a, tstart, tstop, miny, maxy;
  int      tgo_nr, iter_nr, neigh_nr;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  datfile[0]=fitfile[0]=outfile[0]=svgfile[0]=resfile[0]=(char)0;
  dftInit(&dft); fitInit(&fit); resInit(&res);
  tpcProgramName(argv[0], 1, 1, progname, 256);
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    if(strcasecmp(cptr, "W1")==0) {
      weights=1; continue;
    } else if(strcasecmp(cptr, "WF")==0) {
      weights=2; continue;
    } else if(strncasecmp(cptr, "N=", 2)==0) {
      cptr+=2; if(strlen(cptr)>0) {
        is_n_fixed=1; fixed_n=atof_dpi(cptr);
        if(fixed_n>0.0) continue;
      }
    } else if(strncasecmp(cptr, "A=", 2)==0) {
      cptr+=2; if(strlen(cptr)>0) {
        is_a_fixed=1; fixed_a=atof_dpi(cptr);
        if(fixed_a>0.0) continue;
      }
    } else if(strncasecmp(cptr, "Top=", 4)==0) {
      cptr+=4; if(strlen(cptr)>0) {
        is_a_fixed=1; fixed_a=atof_dpi(cptr);
        if(fixed_a>0.0) continue;
      }
    } else if(strncasecmp(cptr, "B=", 2)==0) {
      cptr+=2; if(strlen(cptr)>0) {
        is_b_fixed=1; fixed_b=atof_dpi(cptr);
        parNr=4;
        continue;
      }
    } else if(strncasecmp(cptr, "BOTTOM=", 7)==0) {
      cptr+=7; if(strlen(cptr)>0) {
        is_b_fixed=1; fixed_b=atof_dpi(cptr); parNr=4;
        continue;
      }
    } else if(strcasecmp(cptr, "B")==0) {
      is_b_fixed=0; parNr=4;
      continue;
    } else if(strcasecmp(cptr, "BOTTOM")==0) {
      is_b_fixed=0; parNr=4;
      continue;
    } else if(strcasecmp(cptr, "EC50")==0) {
      parNr=4; type=2801;
      continue;
    } else if(strcasecmp(cptr, "MRL")==0) {
      MRL_check=1; continue;
    } else if(strncasecmp(cptr, "SVG=", 4)==0 && strlen(cptr)>4) {
      strcpy(svgfile, cptr+4); continue;
    } else if(strncasecmp(cptr, "FIT=", 4)==0 && strlen(cptr)>4) {
      strcpy(outfile, cptr+4); continue;
    } else if(strncasecmp(cptr, "RES=", 4)==0 && strlen(cptr)>4) {
      strcpy(resfile, cptr+4); continue;
    }
    fprintf(stderr, "Error: invalid option '%s'\n", argv[ai]);
    return(1);
  } else break;
  
  /* 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 */
  for(; ai<argc; ai++) {
    if(!datfile[0]) {strcpy(datfile, argv[ai]); continue;}
    else if(!fitfile[0]) {strcpy(fitfile, argv[ai]); continue;}
    fprintf(stderr, "Error: invalid argument '%s'\n", argv[ai]);
    return(1);
  }
  if(!fitfile[0]) strcpy(fitfile, "stdout");
  if(type==841 && parNr==4) type=844;
  if(type==841) func=func841;
  else if(type==844) func=func844;
  else func=func2801;
 
  /* Is something missing? */
  if(!datfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  if(is_a_fixed && type!=2801) {
    fprintf(stderr, "Error: option -a can not be used without option -EC50.\n");
    return(1);
  }
  if(is_a_fixed && is_b_fixed && fixed_a<=fixed_b) {
    fprintf(stderr, "Error: invalid constraints for Top and Bottom.\n");
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("datfile := %s\n", datfile);
    printf("fitfile := %s\n", fitfile);
    printf("resfile := %s\n", resfile);
    printf("outfile := %s\n", outfile);
    printf("svgfile := %s\n", svgfile);
    if(is_a_fixed) printf("fixed_a := %g\n", fixed_a);
    if(is_b_fixed) printf("fixed_b := %g\n", fixed_b);
    if(is_n_fixed) printf("fixed_n := %g\n", fixed_n);
    printf("type := %d\n", type);
    printf("weights := %d\n", weights);
    printf("MRL_check := %d\n", MRL_check);
    fflush(stdout);
  }
  if(verbose>4) MATHFUNC_TEST=verbose-4; else MATHFUNC_TEST=0;
 
 
  /*
   *  Read data
   */
  if(verbose>1) printf("reading %s\n", datfile);
  if(dftRead(datfile, &dft)) {
    fprintf(stderr, "Error in reading '%s': %s\n", datfile, dfterrmsg);
    return(2);
  }
  dataNr=dft.frameNr;
 
  /* Sort the samples by time in case data is catenated from several curves */
  (void)dftSortByFrame(&dft);
  if(verbose>30) dftPrint(&dft);
 
  /* Get min and max X and Y */
  ret=dftMinMax(&dft, &tstart, &tstop, &miny, &maxy);
  if(ret!=0) {
    fprintf(stderr, "Error: invalid contents in %s\n", datfile);
    dftEmpty(&dft); return(2);
  }
  if(tstop<=0.0 || maxy<=0.0) {
    fprintf(stderr, "Error: invalid contents in %s\n", datfile);
    dftEmpty(&dft); return(2);
  }
  if(tstart<0.0) fprintf(stderr, "Warning: negative x value(s).\n");
  if(miny<0.0) fprintf(stderr, "Warning: negative y value(s).\n");
 
  /* Check and set weights */
  if(weights==0) {
    if(dft.isweight==0) {
      dft.isweight=1; for(fi=0; fi<dft.frameNr; fi++) dft.w[fi]=1.0;}
  } else if(weights==1) {
    dft.isweight=1; for(fi=0; fi<dft.frameNr; fi++) dft.w[fi]=1.0;
  } else if(weights==2) {
    dftWeightByFreq(&dft);
  }
  if(verbose==10) dftPrint(&dft);
 
 
  /*
   *  Allocate memory for fits
   */
  if(verbose>1) printf("allocating memory for fits.\n");
  ret=fit_allocate_with_dft(&fit, &dft);
  if(ret) {
    fprintf(stderr, "Error: cannot allocate space for fit results (%d).\n", ret);
    dftEmpty(&dft); return(4);
  }
  /* Set necessary information in fit data structure */
  fit.voiNr=dft.voiNr;
  strncpy(fit.datafile, datfile, FILENAME_MAX);
  tpcProgramName(argv[0], 1, 1, fit.program, 1024);
  fit.time=time(NULL);
  for(ri=0; ri<dft.voiNr; ri++) {
    fit.voi[ri].type=type; fit.voi[ri].parNr=parNr;
    fit.voi[ri].start=tstart; fit.voi[ri].end=tstop;
    fit.voi[ri].dataNr=dataNr;
  }
 
 
  /*
   *  Fit the curves
   */
  if(verbose>1) {printf("fitting\n"); fflush(stdout);}
  for(ri=0; ri<dft.voiNr; ri++) {
    if(verbose>0 && dft.voiNr>1) printf("%s\n", dft.voi[ri].name);
 
    /* Set weights for this curve, ie, if NaN, then use 0 weight */
    for(fi=0; fi<dft.frameNr; fi++)
      if(isnan(dft.voi[ri].y[fi])) dft.voi[ri].y3[fi]=0.0;
      else dft.voi[ri].y3[fi]=dft.w[fi];
 
    /* Set data pointers */
    x=dft.x; ymeas=dft.voi[ri].y; yfit=dft.voi[ri].y2; w=dft.voi[ri].y3;
    if(verbose>2) {
      printf("data_weights[%d] := %g", ri+1, dft.w[0]);
      for(fi=1; fi<dft.frameNr; fi++) printf(", %g", w[fi]);
      printf("\n");
    }
 
    /* Set the initial values for parameters */
    fit.voi[ri].wss=3.402823e+38;
 
    /* for defining constraints, find the min and max value of the */
    /* beginning and end of the data set */
    n=dataNr/3; if(n<2) n=2;
    fi=0; while(isnan(dft.voi[ri].y[fi])) fi++;
    min1=max1=dft.voi[ri].y[fi];
    for(fi=0; fi<n; fi++) if(!isnan(dft.voi[ri].y[fi])) {
      if(min1>dft.voi[ri].y[fi]) min1=dft.voi[ri].y[fi];
      if(max1<dft.voi[ri].y[fi]) max1=dft.voi[ri].y[fi];
    }
    fi=dataNr-1; while(isnan(dft.voi[ri].y[fi])) fi--;
    min2=max2=dft.voi[ri].y[fi];
    for(fi=dataNr-n; fi<dataNr; fi++) if(!isnan(dft.voi[ri].y[fi])) {
      if(min2>dft.voi[ri].y[fi]) min2=dft.voi[ri].y[fi];
      if(max2<dft.voi[ri].y[fi]) max2=dft.voi[ri].y[fi];
    }
    if(verbose>3) printf("min1 := %g\nmax1 := %g\nmin2 := %g\nmax2 := %g\n",
      min1, max1, min2, max2);
 
    /* Set parameter constraints */
    if(type==841 || type==844) {
      pmin[0]=(min2>1.0E-08)?min2:1.0E-08;
      pmax[0]=(10.0*max2>pmin[0])?10.0*max2:pmin[0];
      if(is_n_fixed==0) {
        pmin[1]=1.0E-08;  pmax[1]=20.0; //50.0;
      } else {
        pmin[1]=pmax[1]=fixed_n;
      }
      pmin[2]=1.0E-08;  pmax[2]=1.0E20;
      if(is_n_fixed!=0 && fixed_n==1.0) pmax[2]=4.0*dft.x[dataNr-1];
      if(type==844) {
        if(is_b_fixed==0) {
          pmin[3]=0.0; pmax[3]=0.5*(max1+max2);
          //pmin[3]=(min1>1.0e-6)?min1:1.0e-6; pmin[0]-=pmin[3]; pmin[3]*=0.5;
          //max1*=2.0; pmax[3]=(max1>pmin[3])?max1:pmin[3]; //pmax[0]-=pmin[3];
          pmin[0]-=pmax[3];
        } else {
          //pmin[3]=pmax[3]=fixed_b; pmin[0]-=pmin[3]; pmax[0]-=pmin[3];
          pmin[3]=pmax[3]=fixed_b; pmin[0]-=pmax[3]; pmax[0]-=pmin[3];
        }
        if(pmin[0]<0.0) pmin[0]=0.0;
      }
      /* Log10 scaling of parameter 3 */
      pmin[2]=log10(pmin[2]);
      pmax[2]=log10(pmax[2]);
    } else { // type==2801
      /* Top */
      if(is_a_fixed!=0) {
        pmin[0]=pmax[0]=fixed_a; 
      } else {
        a=0.5*(min1+min2);
        pmin[0]=(a>1.0E-10)?a:1.0E-10;
        pmax[0]=10.0*max2;
      }
      /* Bottom */
      if(is_b_fixed!=0) {
        pmin[1]=pmax[1]=fixed_b; 
      } else {
        pmin[1]=0.0; pmax[1]=0.5*(max1+max2);
      }
      /* EC50 */
      pmin[2]=1.0E-08; pmax[2]=4.0*tstop;
      /* HillSlope */
      if(is_n_fixed!=0) {
        pmin[3]=pmax[3]=fixed_n; 
      } else {
        pmin[3]=1.0E-08; pmax[3]=20.0;
      }
    }
    if(verbose>3) {
      fflush(stdout); printf("Constraints for the fit:\n");
      for(pi=0; pi<parNr; pi++)
        printf("  %g - %g\n", pmin[pi], pmax[pi]);
      fflush(stdout);
    }
 
    /* Check that the actual sample number exceeds the number of actually
       fitted parameters */
    for(pi=m=0; pi<parNr; pi++) if(pmin[pi]<pmax[pi]) m++;
    n=getActualSamplenr(&dft, ri);
    if(verbose>2)
      printf("Comparison of nr of samples and params to fit: %d / %d\n", n, m); 
    if((n-1)<m) {
      fprintf(stderr, "Error: too few samples for fitting in %s\n", datfile);
      if(verbose>0) printf("n := %d\nm := %d\n", n, m);
      dftEmpty(&dft); fitEmpty(&fit); return(2);
    }
 
    /* fit */
    if(verbose>2) printf("fitting\n");
    if(type==841) {iter_nr=3; tgo_nr=1000; neigh_nr=15; /*20*/}
    else {iter_nr=3; tgo_nr=2000; neigh_nr=20; /*50*/}
    TGO_LOCAL_INSIDE=0;
    TGO_SQUARED_TRANSF=1;
    ret=tgo(pmin, pmax, func, NULL, parNr, neigh_nr, &fit.voi[ri].wss, 
              fit.voi[ri].p, tgo_nr, iter_nr, verbose-5);
    if(ret) {
      fprintf(stderr, "Error %d in TGO.\n", ret);
      dftEmpty(&dft); fitEmpty(&fit);
      return(4);
    }
 
    /* Warn user about parameters that collide with limits */
    ret=modelCheckLimits(parNr, pmin, pmax, fit.voi[ri].p);
    if(ret!=0 && verbose>0) {
      fprintf(stderr, "Warning: fit collided with the limits.\n");
    }
    /* Force parameter values inside limits, as they are in the simulation */ 
    modelCheckParameters(parNr, pmin, pmax,
                         fit.voi[ri].p, fit.voi[ri].p, &a);
    fit.voi[ri].wss/=a; // remove penalty from reported WSS
 
    /* Print measured and fitted curve */
    func(parNr, fit.voi[ri].p, NULL);
    if(verbose>2) {
      printf("                  Measured  Fitted     Weight\n");
      for(fi=0; fi<dataNr; fi++)
        if(!isnan(ymeas[fi]))
          printf("  %2d: %9.2e   %9.2e %9.2e  %7.1e\n",
                  fi+1, x[fi], ymeas[fi], yfit[fi], w[fi]);
        else
          printf("  %2d: %9.2e   %-9.9s %9.2e  %7.1e\n", fi+1, x[fi],
                  ".", yfit[fi], w[fi]);
      printf(" fitted parameters:");
      for(pi=0; pi<parNr; pi++) printf(" %g", fit.voi[ri].p[pi]); 
      printf("\n");
    }
 
    /* Check the MRL */
    if(verbose>1) {printf("Checking the MRL.\n"); fflush(stdout);}
    m=mrl_between_tacs(yfit, ymeas, dataNr);
    if(m>3 && m>dataNr/3) {
      if(MRL_check!=0) {
        fprintf(stderr, "Error: bad fit.\n");
        fprintf(stderr, "MRL := %d / %d\n", m, dataNr);
        dftEmpty(&dft); fitEmpty(&fit);
        fflush(stderr);
        return(7);
      }
    }
    if(m>2 && m>dataNr/4) {
      fprintf(stderr, "Warning: bad fit.\n");
      fprintf(stderr, "MRL := %d / %d\n", m, dataNr);
      fflush(stderr);
    } else if(verbose>1) {
      printf("MRL test passed.\n"); fflush(stdout);
      if(verbose>2) {
        fprintf(stderr, "MRL := %d / %d\n", m, dataNr); fflush(stderr);}
    }
 
    /* Back-scale p2 */
    if(type==841 || type==844) {
      fit.voi[ri].p[2]=pow(10., fit.voi[ri].p[2]);
    }
 
  } // next curve
 
  /*
   *  Save fit results
   */
  if(verbose>1) printf("saving results in %s\n", fitfile);
  ret=fitWrite(&fit, fitfile); if(ret) {
    fprintf(stderr, "Error (%d) in writing '%s': %s\n", ret, fitfile, fiterrmsg);
    dftEmpty(&dft); fitEmpty(&fit); return(11);
  }
  if(verbose>0) printf("Function parameters written in %s\n", fitfile);
 
  /* Also allocate memory for result file, if needed */
  if(resfile[0]) {
    if(verbose>1) printf("allocating memory for results.\n");
    resInit(&res);
    ret=fitToResult(&fit, &res, temp);
    if(ret!=0) {
      fprintf(stderr, "Error in making results: %s\n", temp);
      dftEmpty(&dft); fitEmpty(&fit); resEmpty(&res); return(4);
    }
    n=0; strcpy(res.parname[n++], "Function");
    if(type==841 || type==844) {
      strcpy(res.parname[n++], "A");
      strcpy(res.parname[n++], "n");
      strcpy(res.parname[n++], "K");
    }
    if(type==844) {
      strcpy(res.parname[n++], "B");
    }
    if(type==2801) {
      strcpy(res.parname[n++], "Top");
      strcpy(res.parname[n++], "Bottom");
      strcpy(res.parname[n++], "EC50");
      strcpy(res.parname[n++], "HillSlope");
    }
    strcpy(res.parname[n++], "WSS");
    strcpy(res.program, progname);
    sprintf(res.datarange, "%g - %g", tstart, tstop);
    strcpy(res.studynr, dft.studynr);
    if(verbose>0) resPrint(&res);
    if(verbose>1) printf("saving results in %s\n", resfile);
    if(resWrite(&res, resfile, verbose-3)) {
      fprintf(stderr, "Error in writing '%s': %s\n", resfile, reserrmsg);
      dftEmpty(&dft); fitEmpty(&fit); resEmpty(&res); return(11);
    }
    resEmpty(&res);
    if(verbose>1) printf("Function parameters written in %s\n", resfile);
  }
 
 
  /*
   *  Plot fitted curves
   */
  if(svgfile[0]) {
    if(verbose>1) printf("creating SVG plot\n");
    DFT adft;
    char tmp[FILENAME_MAX];
 
    if(verbose>1)
      printf("calculating fitted curve at automatically generated sample times\n");
    dftInit(&adft);
    ret=dftAutointerpolate(&dft, &adft, 1.05*tstop, verbose-10);
    if(ret) {
      fprintf(stderr, "Error %d in memory allocation for fitted curves.\n", ret);
      dftEmpty(&dft); fitEmpty(&fit); dftEmpty(&adft);
      return(13);
    }
    //for(fi=0; fi<adft.frameNr; fi++) adft.w[fi]=1.0;
    for(ri=0, ret=0; ri<adft.voiNr; ri++) {
      for(fi=0; fi<adft.frameNr; fi++) {
        ret=fitEval(&fit.voi[ri], adft.x[fi], &a); if(ret!=0) break;
        adft.voi[ri].y[fi]=a;
      }
      if(ret!=0) break;
    }
    if(ret!=0) {
      fprintf(stderr, "Error: cannot calculate fitted TAC for '%s'.\n", svgfile);
      dftEmpty(&dft); dftEmpty(&adft); fitEmpty(&fit);
      return(14);
    }
 
    /* Save SVG plot of fitted and original data */
    if(verbose>1) printf("writing %s\n", svgfile);
    tpcProgramName(argv[0], 0, 0, tmp, 128); strcat(tmp, " ");
    if(strlen(adft.studynr)>1) strcat(tmp, adft.studynr);
    ret=plot_fitrange_svg(&dft, &adft, tmp, 0.0, nan(""), 0.0, nan(""),
                          svgfile, verbose-10);
    if(ret) {
      fprintf(stderr, "Error (%d) in writing '%s'.\n", ret, svgfile);
      dftEmpty(&adft); dftEmpty(&dft); fitEmpty(&fit);
      return(30+ret);
    }
    if(verbose>0) printf("Plots written in %s\n", svgfile);
 
    dftEmpty(&adft);
  }
 
 
  /*
   *  Save fitted TACs (original data is lost here, so do this in the end)
   */
  if(outfile[0]) {
    if(verbose>1) printf("saving fitted curves in %s\n", outfile);
 
    for(ri=0, ret=0; ri<dft.voiNr; ri++) {
      for(fi=0; fi<dft.frameNr; fi++) {
        ret=fitEval(&fit.voi[ri], dft.x[fi], &a); if(ret!=0) break;
        dft.voi[ri].y[fi]=a;
      }
      if(ret!=0) break;
    }
    if(ret!=0) {
      fprintf(stderr, "Error: cannot calculate fitted TACs for '%s'.\n", outfile);
      dftEmpty(&dft); fitEmpty(&fit);
      return(21);
    }
    tpcProgramName(argv[0], 1, 0, temp, 128);
    sprintf(dft.comments, "# program := %s\n", temp);
    dft.isweight=0;
    if(dftWrite(&dft, outfile)) {
      fprintf(stderr, "Error in writing '%s': %s\n", outfile, dfterrmsg);
      dftEmpty(&dft); fitEmpty(&fit);
      return(22);
    }
    if(verbose>0) printf("Fitted TACs written in %s\n", outfile);
  }
 
  /* Free memory */
  dftEmpty(&dft); fitEmpty(&fit);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************
 *
 *  Functions to be minimized
 *
 *****************************************************************************/
double func841(int parNr, double *p, void *fdata)
{
  int i;
  double a, b, c, q, v, wss;
  double pa[MAX_PARAMETERS], penalty=1.0;
 
  /* Check parameters against the constraints */
  (void)modelCheckParameters(parNr, pmin, pmax, p, pa, &penalty);
  if(fdata) {}
 
  /* Get parameters */
  a=pa[0]; b=pa[1]; c=pow(10., pa[2]);
  /* Calculate the curve values and weighted SS */
  for(i=0, wss=0.0; i<dataNr; i++) {
    if(x[i]>=0.0) q=pow(x[i], b); else q=1.0;
    yfit[i]= a*q/(q+c); if(isnan(ymeas[i])) continue;
    v=yfit[i]-ymeas[i]; wss+=w[i]*v*v;
  }
  // printf("%g %g %g -> %g\n", a, b, c, wss);
  wss*=penalty;
 
  return(wss);
}
 
double func844(int parNr, double *p, void *fdata)
{
  int i;
  double a, b, c, d, q, v, wss;
  double pa[MAX_PARAMETERS], penalty=1.0;
  if(fdata) {}
 
  /* Check parameters against the constraints */
  (void)modelCheckParameters(parNr, pmin, pmax, p, pa, &penalty);
 
  /* Get parameters */
  a=pa[0]; b=pa[1]; c=pow(10., pa[2]); d=pa[3];
  /* Calculate the curve values and weighted SS */
  for(i=0, wss=0.0; i<dataNr; i++) {
    if(ymeas[i]<min1) continue;
    if(x[i]>=0.0) q=pow(x[i], b); else q=1.0;
    yfit[i]= a*q/(q+c) + d; if(isnan(ymeas[i])) continue;
    v=yfit[i]-ymeas[i]; wss+=w[i]*v*v;
  }
  // printf("%g %g %g %g -> %g\n", a, b, c, d, wss);
  wss*=penalty;
 
  return(wss);
}
 
double func2801(int parNr, double *p, void *fdata)
{
  int i;
  double top, bottom, ec50, hillslope, v, wss;
  double pa[MAX_PARAMETERS], penalty=1.0;
  if(fdata) {}
 
  /* Check parameters against the constraints */
  (void)modelCheckParameters(parNr, pmin, pmax, p, pa, &penalty);
 
  /* Get parameters */
  top=pa[0]; bottom=pa[1]; ec50=pa[2]; hillslope=pa[3];
  /* Calculate the curve values and weighted SS */
  for(i=0, wss=0.0; i<dataNr; i++) {
    if(x[i]<0.0) {yfit[i]=bottom; continue;}
    if(x[i]<1.0E-12) yfit[i]=bottom;
    else yfit[i]=bottom+(top-bottom)/(1.0+pow(ec50/x[i],hillslope));
    if(isnan(ymeas[i])) continue;
    v=yfit[i]-ymeas[i]; wss+=w[i]*v*v;
  }
  wss*=penalty;
 
  return(wss);
}
/*****************************************************************************/
 
/*****************************************************************************/