fit_fexp.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"
/*****************************************************************************/
double *x, *ymeas, *yfit, *w;
int fitdataNr=0, parNr=3;
double pmin[MAX_PARAMETERS], pmax[MAX_PARAMETERS];
double lastWSS;
double func351(int parNr, double *p, void*);
double funcMono(int parNr, double *p, void*);
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Fits the following exponential function to plasma parent (unchanged) tracer",
  "fraction curve, where fractions are between 0 and 1 (1,2).",
  " ",
  "                             -B*x  -C*x                               ",
  "        f(x) = 1 - A * (2 - e    -e    )   , where 0<A<=1 , B>0 , C>0 ",
  " ",
  "Optionally a monoexponential function, where fractions are between A and B,",
  "can be applied instead:",
  "                        -C*x                               ",
  "        f(x) = (A-B) * e      + B          , where 0<A<=1 , B>=0 , C>0 ",
  " ",
  "Usage: @P [Options] fractionfile [fitfile]",
  " ",
  "Options:",
  " -mono",
  "     Use monoexponential function.",
  " -a=<value> or -a=zero",
  "     Parent fraction at time zero (A) is fixed to given value, or to",
  "     the value of zero sample, if available.",
  "     Applicable only with monoexponential function.",
  " -WC",
  "     The last data column contains sample weights.",
  " -W1",
  "     All weights are set to 1.0 (no weighting)",
  " -WF",
  "     Less frequent samples are given more weight.",
  " -ND",
  "     Some fractions are known to exceed 1, not divided by 100.",
  " -MRL",
  "     Error is returned if MRL check is not passed.",
  " -svg=<Filename>",
  "     Fitted and measured TACs are plotted in specified SVG file.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Fraction datafile must contain two columns, sample times (min) and fractions",
  "of parent tracer. Weights can be specified as specified in DFT format (4);",
  "any additional columns are ignored.",
  "Lines that start with a '#' are not read from the datafile.",
  "Program writes the fit start and end times, nr of points, WSS/n,",
  "and parameters of the fitted function to the FIT file (5).",
  " ",
  "References:",
  "1. Fitting the fractions of parent tracer in plasma.",
  "   https://www.turkupetcentre.net/petanalysis/input_parent_fitting.html",
  "2. Kropholler MA, Boellaard R, Schuitemaker A, van Berckel BNM, Luurtsema G,",
  "   Windhorst AD, Lammertsma AA. Development of a tracer kinetic plasma input",
  "   model for (R)-[11C]PK11195 brain studies. J Cereb Blood Flow Metab. 2005;",
  "   25(7): 842-851.",
  "4. https://www.turkupetcentre.net/petanalysis/format_tpc_dft.html",
  "5. https://www.turkupetcentre.net/petanalysis/format_tpc_fit.html",
  " ",
  "See also: fit2dat, metabcor, avgfract, fit_ppf, fit_dexp, tac2svg",
  " ",
  "Keywords: input, plasma, modelling, simulation, metabolite correction",
  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=351, ret, m, n;
  int        noDivide=0;
  int        MRL_check=0; // 0=off, 1=on
  double     fixed_a=nan("");
  int        fixed0=0; // fix fraction at time zero to measured value at time zero, if available
  char      *cptr; 
  char       dfile[FILENAME_MAX], ffile[FILENAME_MAX], svgfile[FILENAME_MAX];
  DFT        dft;
  FIT        fit;
  double     a, tstart, tstop, d, maxfract, aic;
  // Weights: 0=left as it is in datafile, 1=set to 1, 2=sampling frequency
  int        weighting=0; 
  int        last_col_is_weight=0;
 
#ifdef MINGW
  // Use Unix/Linux default of two-digit exponents in MinGW on Windows
  _set_output_format(_TWO_DIGIT_EXPONENT);
#endif
 
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  dftInit(&dft); fitInit(&fit);
  dfile[0]=ffile[0]=svgfile[0]=(char)0;
  /* 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, "ND")==0) {
      noDivide=1; continue;
    } else if(strcasecmp(cptr, "MONO")==0) {
      type=302; parNr=4; continue;
    } else if(strcasecmp(cptr, "A=ZERO")==0 || strcasecmp(cptr, "A=0")==0) {
      fixed0=1; continue;
    } else if(strncasecmp(cptr, "A=", 2)==0) {
      fixed_a=atof_dpi(cptr+2); if(fixed_a>0.0) continue;
    } else if(strcasecmp(cptr, "W1")==0) {
      weighting=1; continue;
    } else if(strcasecmp(cptr, "WF")==0) {
      weighting=2; continue;
    } else if(strcasecmp(cptr, "WC")==0) {
      last_col_is_weight=1; continue;
    } else if(strcasecmp(cptr, "MRL")==0) {
      MRL_check=1; continue;
    } else if(strncasecmp(cptr, "SVG=", 4)==0 && strlen(cptr)>4) {
      strlcpy(svgfile, cptr+4, FILENAME_MAX); continue;
    }
    /* we should never get this far */
    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 */
  if(ai<argc) strlcpy(dfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(ffile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {fprintf(stderr, "Error: too many arguments: '%s'.\n", argv[ai]); return(1);}
 
  /* Is something missing? */
  if(!dfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  if(!ffile[0]) strcpy(ffile, "stdout");
  if(verbose>6) MATHFUNC_TEST=verbose-6; else MATHFUNC_TEST=0;
  if(verbose>6) CSV_TEST=verbose-6; else CSV_TEST=0;
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("dfile := %s\n", dfile);
    printf("ffile := %s\n", ffile);
    printf("svgfile := %s\n", svgfile);
    printf("noDivide := %d\n", noDivide);
    printf("weighting := %d\n", weighting);
    printf("last_col_is_weight := %d\n", last_col_is_weight);
    printf("fixed0 := %d\n", fixed0);
    printf("MRL_check := %d\n", MRL_check);
    printf("type := %d\n", type);
    printf("parNr := %d\n", parNr);
  }
 
 
  /*
   *  Read data
   */
  if(verbose>1) printf("reading %s\n", dfile);
  if(dftRead(dfile, &dft)) {
    fprintf(stderr, "Error in reading '%s': %s\n", dfile, dfterrmsg);
    return(2);
  }
  if(dft.frameNr<2) { // More precise testing later
    fprintf(stderr, "Error: too few samples for fitting in %s\n", dfile);
    dftEmpty(&dft); return(2);
  }
  fitdataNr=dft.frameNr;
 
  /* If required, set the last column as weights */
  if(last_col_is_weight!=0) {
    if(verbose>1) printf("reading weights from the last column.\n");
    if(dft.voiNr<2) {
      fprintf(stderr, "Error: no column for weights in %s\n", dfile);
      dftEmpty(&dft); return(2);
    }
    for(fi=0; fi<dft.frameNr; fi++) {
      dft.w[fi]=dft.voi[dft.voiNr-1].y[fi];
      if(isnan(dft.w[fi])) dft.w[fi]=0.0;
    }
    dft.isweight=1; dft.voiNr--;
  }
 
  if(verbose>1) {
    if(strlen(dft.studynr)>0 && strcmp(dft.studynr, ".")!=0)
      printf("study_number := %s\n", dft.studynr);
    printf("tacNr := %d\n", dft.voiNr);
  }
 
 
  /* Ignore extra columns */
  if(dft.voiNr>1) {
    fprintf(stderr, "Warning: extra columns in %s are ignored.\n", dfile);
    dft.voiNr=1;
  }
 
  /* Sort the samples by time in case data is catenated from several curves */
  (void)dftSortByFrame(&dft);
 
  /* Make sure that times are in minutes */
  if(dft.timeunit==TUNIT_SEC) {
    if(verbose>1) printf("Sample times are converted to min\n");
    dftSec2min(&dft);
  }
  if(dft.timeunit==TUNIT_UNKNOWN) {
    dft.timeunit=TUNIT_MIN;
  }
 
  /* Set the names for TAC */
  if(strlen(dft.voi[0].name)==0) {
    strcpy(dft.voi[0].name, "Parent");
    strcpy(dft.voi[0].voiname, dft.voi[0].name);
  }
 
  /* Get start and end times */
  ret=dftMinMax(&dft, &tstart, &tstop, NULL, &maxfract);
  if(ret!=0) {
    fprintf(stderr, "Error: invalid contents in %s\n", dfile);
    dftEmpty(&dft); return(2);
  }
  /* Check that x values are >=0 and that there is range in x values */
  if(tstart<0.0 || !(tstop>tstart)) {
    fprintf(stderr, "Error: invalid sample times.\n");
    dftEmpty(&dft); return(2);
  }
  /* Check that y values are <=1 */
  if(maxfract>100.0 || maxfract<0.001) {
    fprintf(stderr, "Error: invalid fraction values.\n");
    dftEmpty(&dft); return(2);
  }
  if(noDivide==0 && maxfract>1.0) {
    fprintf(stderr, "Warning: converting percentages to fractions.\n");
    for(fi=0; fi<dft.frameNr; fi++) for(ri=0; ri<dft.voiNr; ri++)
      if(!isnan(dft.voi[ri].y[fi])) dft.voi[ri].y[fi]/=100.0;
    /* make sure that units are not percentages any more */
    dftUnitToDFT(&dft, CUNIT_UNITLESS);
  }
  if(verbose>9) dftPrint(&dft);
 
  /* Check and set weights */
  if(dft.isweight==0 || weighting==1) {
    dft.isweight=1;
    for(fi=0; fi<dft.frameNr; fi++) dft.w[fi]=1.0;
  }
  if(weighting==2) { // weighting by sample frequency
    dftWeightByFreq(&dft); // must be sorted by time before this
  }
  /* Set sample weights to zero, if parent fraction is NaN */
  for(fi=0; fi<dft.frameNr; fi++) if(isnan(dft.voi[0].y[fi])) dft.w[fi]=0;
  if(verbose>1) {
    printf("data_weights := %g", dft.w[0]);
    for(fi=1; fi<dft.frameNr; fi++) printf(", %g", dft.w[fi]);
    printf("\n");
  }
 
  /* If there is zero sample, then get its value */
  double y0=nan("");
  if(fabs(tstart)<1.0E-08) {
    /* There might be more than one zero samples, thus get the mean */
    int n=0; y0=0.0;
    for(int i=0; i<fitdataNr; i++) if(fabs(dft.x[i])<1.0E-08) {y0+=dft.voi[0].y[i]; n++;}
    if(n>1) y0/=(double)n;
    if(y0>0.0 && y0<=1.0) {
      if(verbose>1) {printf("initial_fraction := %g\n", y0); fflush(stdout);}
    } else {
      y0=nan("");
    }
  }
 
 
  /*
   *  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); fitEmpty(&fit); return(4);
  }
  /* Set necessary information in fit data structure */
  fit.voiNr=dft.voiNr;
  strlcpy(fit.datafile, dfile, FILENAME_MAX);
  tpcProgramName(argv[0], 1, 1, fit.program, 256);
  fit.time=time(NULL);
  for(fi=n=0; fi<fitdataNr; fi++)
    if(dft.x[fi]>=0.0 && !isnan(dft.voi[0].y[fi])) n++;
  for(ri=0; ri<fit.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=n;
  }
 
 
  /*
   *  Fit
   */
  if(verbose>1) printf("fitting\n");
  ri=0;
  /* Set data pointers */
  x=dft.x; ymeas=dft.voi[ri].y; yfit=dft.voi[ri].y2; w=dft.w;
  /* Set the initial values for parameters */
  fit.voi[ri].wss=3.402823e+38;
 
  /* Set parameter constraints */
  if(type==302) {
    pmin[0]=0.0;  pmax[0]=1.0;
    pmin[1]=0.0;  pmax[1]=1.0;
    pmin[2]=0.0;  pmax[2]=1.5;
    if(!isnan(fixed_a)) {pmin[0]=pmax[0]=fixed_a; pmax[1]=fixed_a;}
    else if(fixed0 && !isnan(y0)) {pmin[0]=pmax[0]=y0; pmax[1]=y0;}
  } else {
    pmin[0]=0.0;  pmax[0]=0.5;
    pmin[1]=0.0;  pmax[1]=1.0;
    pmin[2]=0.0;  pmax[2]=0.5;
  }
  if(verbose>2) {
    printf("Constraints for the 1st fit:\n");
    for(pi=0; pi<parNr; pi++)
      printf("  %g - %g\n", pmin[pi], pmax[pi]);
  }
 
 
  /*
   *  Check that sample number (not including NaNs) is at least one more
   *  than the number of actually fitted parameters.
   */
  for(fi=n=0; fi<fitdataNr; fi++) {
    if(dft.x[fi]<0.0) continue;
    if(dft.w[fi]>0.0) n++;
  }
  for(pi=m=0; pi<parNr; pi++) if(pmin[pi]<pmax[pi]) m++;
  if(verbose>1) 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", dfile);
    if(verbose>0) printf("n := %d\nm := %d\n", n, m);
    dftEmpty(&dft);  fitEmpty(&fit); return(2);
  }
 
 
 
  /* fit */
  TGO_LOCAL_INSIDE=1;
  TGO_LOCAL_OPT=0; // 0=Powell-Brent, 1=Bobyqa
  TGO_SQUARED_TRANSF=1; 
  int neighNr=5, iterNr=1, sampleNr=1000;
  if(type==351) 
    ret=tgo(pmin, pmax, func351, NULL, parNr, neighNr, &fit.voi[ri].wss, 
            fit.voi[ri].p, sampleNr, iterNr, verbose-8);
  else
    ret=tgo(pmin, pmax, funcMono, NULL, parNr, neighNr, &fit.voi[ri].wss, 
            fit.voi[ri].p, sampleNr, iterNr, verbose-8);
  if(ret) {
    fprintf(stderr, "Error %d in TGO.\n", ret);
    dftEmpty(&dft); fitEmpty(&fit);
    return(4);
  }
  if(verbose>4) {
    printf("Results from tgo():\n");
    for(pi=0; pi<parNr; pi++)
      printf("  parameter[%d] := %g\n", pi, fit.voi[ri].p[pi]);
    printf("WSS := %g\n", fit.voi[ri].wss);
  }
  fit.voi[ri].wss=lastWSS; // remove penalty from reported WSS (do not use a)
  if(verbose>5) printf("lastWSS := %g\n", lastWSS);
  /* Correct fitted parameters to match constraints like inside the function */
  (void)modelCheckParameters(parNr, pmin, pmax, fit.voi[ri].p, fit.voi[ri].p, NULL);
 
  /* Further processing with bi-exponential */
  if(type==351) {
 
    /* Sort b and c */
    if(fit.voi[ri].p[1]<fit.voi[ri].p[2]) {
      a=fit.voi[ri].p[1]; fit.voi[ri].p[1]=fit.voi[ri].p[2]; fit.voi[ri].p[2]=a;
      a=pmin[1]; pmin[1]=pmin[2]; pmin[2]=a;
      a=pmax[1]; pmax[1]=pmax[2]; pmax[2]=a;
    }
 
    /* If b and c are about the same, fit again only a and b, with c=0 */
    d=fabs(fit.voi[ri].p[1]-fit.voi[ri].p[2]);
    a=0.5*(fit.voi[ri].p[1]+fit.voi[ri].p[2]);
    if(verbose>1) printf("|b-c|=%g mean(b, c)=%g\n", d, a);
    if(a>0.0 && d/a<0.01) {
      if(verbose>1) printf("Fitting again with only one exponential...\n");
      pmin[0]=0.0;  pmax[0]=1.0;
      pmin[1]=0.0;  pmax[1]=1.5;
      pmin[2]=0.0;  pmax[2]=0.0;
      if(verbose>2) {
        printf("Constraints for the 2nd fit:\n");
        for(pi=0; pi<parNr; pi++)
          printf("  %g - %g\n", pmin[pi], pmax[pi]);
      }
      neighNr=4, iterNr=1, sampleNr=1000;
      ret=tgo(pmin, pmax, func351, NULL, parNr, neighNr, &fit.voi[ri].wss, 
              fit.voi[ri].p, sampleNr, iterNr, verbose-8);
      if(ret) {
        fprintf(stderr, "Error %d in TGO.\n", ret);
        dftEmpty(&dft); fitEmpty(&fit);
        return(4);
      }
      fit.voi[ri].wss=lastWSS; // non-penalized WSS
      /* Correct fitted parameters to match constraints like inside the function */
      (void)modelCheckParameters(parNr, pmin, pmax, fit.voi[ri].p, fit.voi[ri].p, NULL);
    }
    /* Warn user about parameters that collide with limits */
    if(verbose>1) {
      ret=modelCheckLimits(parNr, pmin, pmax, fit.voi[ri].p);
      if(ret==0) {if(verbose>2) fprintf(stdout, "ok\n");}
      else fprintf(stderr, "warning, fit collided with the limits.\n");
    }
    if(verbose>1) fprintf(stdout, "WSS := %g\n", fit.voi[ri].wss);
  }
 
 
  /* Further processing with monoexponential */
  if(type==302) {
    // Function has four parameters, and order is not the same
    // f(x)=(A-B)*exp(-C*x)+B
    // f(x)=A*exp(B*x)+C*exp(D*x)
    double a=fit.voi[ri].p[0]-fit.voi[ri].p[1];
    double b=-fit.voi[ri].p[2];
    double c=fit.voi[ri].p[1];
    double d=0.0;
    fit.voi[ri].p[0]=a;
    fit.voi[ri].p[1]=b;
    fit.voi[ri].p[2]=c;
    fit.voi[ri].p[3]=d;
  }
 
  /* Print measured and fitted curve */
  if(verbose>3) {
    printf("     Measured  Fitted:\n");
    for(fi=0; fi<fitdataNr; fi++)
      if(!isnan(ymeas[fi]))
        printf("  %2d: %9.2e   %9.2e %9.2e\n", fi+1, x[fi], ymeas[fi], yfit[fi]);
      else
        printf("  %2d: %9.2e   %-9.9s %9.2e\n", fi+1, x[fi], ".", yfit[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");
  fi=mrl_between_tacs(yfit, ymeas, fitdataNr);
  if(fi>3 && fi>fitdataNr/3) {
    if(MRL_check!=0) {
      fprintf(stderr, "Error: bad fit.\n");
      fprintf(stderr, "MRL := %d / %d\n", fi, fitdataNr);
      dftEmpty(&dft); fitEmpty(&fit);
      return(7);
    }
  }
  if(fi>2 && fi>fitdataNr/4) {
    fprintf(stderr, "Warning: bad fit.\n");
    fprintf(stderr, "MRL := %d / %d\n", fi, fitdataNr);
  } else if(verbose>1) {
    printf("MRL test passed.\n");
    if(verbose>2) fprintf(stderr, "MRL := %d / %d\n", fi, fitdataNr);
  }
 
 
  /* Calculate AIC */
  if(verbose>1) {
    int m;
    for(fi=0, m=0; fi<fitdataNr; fi++) if(w[fi]>0.0) m++;
    for(pi=0, n=0; pi<parNr; pi++) if(pmax[pi]>pmin[pi]) n++;
    printf("nr_of_fitted_parameters := %d\n", n);
    aic=aicSS(fit.voi[ri].wss, m, n);
    printf("AIC := %g\n", aic);
  }
 
 
  /*
   *  Save fit results
   */
  if(verbose>1 && strcmp(ffile, "stdout")!=0) printf("saving results\n");
  ret=fitWrite(&fit, ffile); if(ret) {
    fprintf(stderr, "Error (%d) in writing '%s': %s\n", ret, ffile, fiterrmsg);
    dftEmpty(&dft); fitEmpty(&fit); return(11);
  }
  if(strcmp(ffile, "stdout")!=0 && verbose>0)
    printf("Function parameters written in %s\n", ffile);
 
 
  /*
   *  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.10*dft.x[dft.frameNr-1], verbose-7);
    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 curve 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, 64); strcat(tmp, " ");
    if(strlen(adft.studynr)>1) strcat(tmp, adft.studynr);
    ret=plot_fitrange_svg(&dft, &adft, tmp, 0.0, nan(""), 0.0, 1.0, svgfile, verbose-8);
    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);
  }
 
 
  /* Free memory */
  dftEmpty(&dft); fitEmpty(&fit);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************
 *
 *  Functions to be minimized
 *
 *****************************************************************************/
double func351(int parNr, double *p, void *fdata)
{
  int i;
  double pa[MAX_PARAMETERS], a, b, c, v, wss, penalty=1.0;
 
  /* Check parameters against the constraints */
  if(0) {
    for(i=0; i<parNr; i++) printf(" p[%d]=%g", i, p[i]);
    printf("\n");
  }
  modelCheckParameters(parNr, pmin, pmax, p, pa, &penalty);
  if(fdata) {}
  a=pa[0]; b=pa[1]; c=pa[2];
  /* Calculate the curve values and weighted SS */
  for(i=0, wss=0.0; i<fitdataNr; i++) if(w[i]>0.0) {
    yfit[i]= 1.0 - a*(2.0-exp(-b*x[i])-exp(-c*x[i]));
    if(yfit[i]<0.0 || yfit[i]>1.0) penalty*=2.0;
    v=yfit[i]-ymeas[i]; wss+=w[i]*v*v;
  }
  lastWSS=wss;
  wss*=penalty;
  if(0) printf("%g %g %g -> %g\n", a, b, c, wss);
  return(wss);
}
 
 
 
double funcMono(int parNr, double *p, void *fdata)
{
  int i;
  double pa[MAX_PARAMETERS], a, b, c, v, wss, penalty=1.0;
 
  /* Check parameters against the constraints */
  if(0) {
    for(i=0; i<parNr; i++) printf(" p[%d]=%g", i, p[i]);
    printf("\n");
  }
  modelCheckParameters(parNr, pmin, pmax, p, pa, &penalty);
  if(fdata) {}
  a=pa[0]; b=pa[1]; c=pa[2];
  /* Calculate the curve values and weighted SS */
  for(i=0, wss=0.0; i<fitdataNr; i++) if(w[i]>0.0) {
    yfit[i]= (a-b)*exp(-c*x[i]) + b;
    if(yfit[i]<0.0 || yfit[i]>1.0) penalty*=2.0;
    v=yfit[i]-ymeas[i]; wss+=w[i]*v*v;
  }
  lastWSS=wss;
  wss*=penalty;
  if(0) printf("%g %g %g -> %g\n", a, b, c, wss);
  return(wss);
}
/*****************************************************************************/
 
/*****************************************************************************/