fit_rrtm.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <unistd.h>
#include <string.h>
#include <time.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcmodel.h"
#include "libtpccurveio.h"
#include "libtpcsvg.h"
#include "libtpcmodext.h"
/*****************************************************************************/
 
/*****************************************************************************/
const int parNr=3;
int fitframeNr=0;
double *t, *cr, *ct, *tis, *w; /* These are pointers, not allocated */
double pmin[MAX_PARAMS], pmax[MAX_PARAMS];
double wss_wo_penalty=0.0;
/* Local functions */
double rrtmFunc(int parNr, double *p, void*);
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "NLLSQ estimation of R1 (=K1/K1'), k2, and k3 applying the reduced reference",
  "tissue compartment model, RRTM (1). This model is based on the reference",
  "tissue compartment model (2, 3), but here it is assumed that the binding or",
  "metabolism of tracer is irreversible (k4=0) during the PET scanning.",
  " ",
  "Usage: @P [Options] ttacfile reference endtime resultfile",
  " ",
  "TTAC file can be in DFT or PMOD format. Sample times must be in minutes.",
  "If TTAC file contains weights, those are used in the NLLSQ fitting.",
  "Reference region TAC can be given separate TAC file or as the name or number",
  "of the reference region in TTAC file.",
  " ",
  "Options:",
  " -lim=<filename>",
  "     Specify the constraints for model parameters;",
  "     This file with default values can be created by giving this option",
  "     as the only command-line argument to this program.",
  " -SD[=<y|N>]",
  "     Standard deviations are calculated and saved in results (y), or",
  "     not calculated (n).",
  " -CL[=<y|N>]",
  "     95% Confidence limits are calculated and saved in results (y), or",
  "     not calculated (n).",
  " -w1",
  "     All weights are set to 1.0 (no weighting); by default, weights in",
  "     TTAC file are used, if available.",
  " -wf",
  "     Weight by sampling interval.",
  " -fit=<Filename>",
  "     Fitted regional TACs are written in file.",
  " -svg=<Filename>",
  "     Fitted and measured TACs are plotted in specified SVG file.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  " ",
  "Values of R1, k2, and k3 are written in the specified result file.",
  "Fitted curves are written in DFT format, if file name is given.",
  " ",
  "Example 1: file a789.dft contains regions-of-interest and reference region,",
  "with name 'cereb'. The whole time range is used in the fit.",
  "Fitted TACs are plotted in SVG format.",
  "     @P -svg=a789rrtm.svg a789.dft cereb 999 a789.res",
  " ",
  "Example 2: Reference region TAC is in a separate file, a789ref.dft.",
  "Standard deviations and confidence limits are also estimated.",
  "TAC data from injection to 60 min is used in the fitting.",
  "     @P -SD=y -CL=y a789.dft a789ref.dft 60 a789.res",
  " ",
  "Example 3a: Create a file containing default parameter limits:",
  "     @P -lim=rrtm.lim",
  " ",
  "Example 3b: Apply user-defined parameter constraints specified in rrtm.lim:",
  "     @P -lim=rrtm.lim a789.dft cereb 999 a789.res",
 
  "References:",
  "1. Oikonen V. Model equations for reference tissue compartmental models.",
  "   https://www.turkupetcentre.net/reports/tpcmod0002.pdf",
  "2. Cunningham VJ, Hume SP, Price GR, Ahier RG, Cremer JE, Jones AKP.",
  "   Compartmental analysis of diprenorphine binding to opiate receptors",
  "   in the rat in vivo and its comparison with equilibrium data in vitro.",
  "   J Cereb Blood Flow Metab 1991;11:1-9.",
  "3. Lammertsma AA, Hume SP. Simplified reference tissue model for PET",
  "   receptor studies. NeuroImage 1996;4:153-158.",
  " ",
  "See also: patlak, fitk3, tacweigh, rescoll, sim_rtcm, fit_trtm",
  " ",
  "Keywords: TAC, modelling, irreversible uptake, RTCM, reference input",
  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    rtacfile[FILENAME_MAX], ttacfile[FILENAME_MAX], resfile[FILENAME_MAX],
          fitfile[FILENAME_MAX], svgfile[FILENAME_MAX], limfile[FILENAME_MAX];
  char   *cptr, tmp[256];
  double  fitdur=nan("");
  int     weights=0; // 0=default, 1=no weighting, 2=frequency
  int     doBootstrap=0, doSD=0, doCL=0;
  double *sd, *cl1, *cl2;
  double  def_pmin[MAX_PARAMETERS], def_pmax[MAX_PARAMETERS];
  int     ret, n;
 
 
  /* Set parameter initial values and constraints */
  /* R1  */ def_pmin[0]=0.001;     def_pmax[0]=2.0;
  /* k2  */ def_pmin[1]=0.000001;  def_pmax[1]=0.5;
  /* k3  */ def_pmin[2]=0.0;       def_pmax[2]=1.0;
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  rtacfile[0]=ttacfile[0]=resfile[0]=fitfile[0]=svgfile[0]=limfile[0]=(char)0;
  /* Get options first, because it affects what arguments are read */
  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(strncasecmp(cptr, "CL", 2)==0) {
      if(strlen(cptr)==2) {doCL=1; continue;}
      cptr+=2; if(*cptr=='=') {
        cptr++;
        if(*cptr=='Y' || *cptr=='y') {doCL=1; continue;}
        if(*cptr=='N' || *cptr=='n') {doCL=0; continue;}
      }
    } else if(strncasecmp(cptr, "SD", 2)==0) {
      if(strlen(cptr)==2) {doSD=1; continue;}
      cptr+=2; if(*cptr=='=') {
        cptr++;
        if(*cptr=='Y' || *cptr=='y') {doSD=1; continue;}
        if(*cptr=='N' || *cptr=='n') {doSD=0; continue;}
      }
    } else if(strncasecmp(cptr, "LIM=", 4)==0 && strlen(cptr)>4) {
      strlcpy(limfile, cptr+4, FILENAME_MAX); continue;
    } else if(strcasecmp(cptr, "LIM")==0) {
      strcpy(limfile, "stdout"); continue;
    } else if(strcasecmp(cptr, "W1")==0) {
      weights=1; continue;
    } else if(strcasecmp(cptr, "WF")==0) {
      weights=2; continue;
    } else if(strncasecmp(cptr, "FIT=", 4)==0) {
      strlcpy(fitfile, cptr+4, FILENAME_MAX); if(strlen(fitfile)>0) continue;
    } else if(strncasecmp(cptr, "SVG=", 4)==0) {
      strlcpy(svgfile, cptr+4, FILENAME_MAX); if(strlen(svgfile)>0) 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 */
  if(ai<argc) strlcpy(ttacfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(rtacfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {
    if(atof_with_check(argv[ai], &fitdur)!=0 || fitdur<0.0) {
      fprintf(stderr, "Error: invalid fit time: '%s'.\n", argv[ai]);
      return(1);
    }
    if(fitdur==0) fitdur=1.0E+10;
    ai++;
  }
  if(ai<argc) strlcpy(resfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
  if(doSD || doCL) doBootstrap=1; else doBootstrap=0;
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("ttacfile := %s\n", ttacfile);
    printf("reference := %s\n", rtacfile);
    printf("resfile := %s\n", resfile);
    printf("fitfile := %s\n", fitfile);
    printf("svgfile := %s\n", svgfile);
    printf("limfile := %s\n", limfile);
    printf("required_fittime := %g min\n", fitdur);
    printf("weights := %d\n", weights);
    printf("doBootstrap := %d\n", doBootstrap);
    printf("doSD := %d\n", doSD);
    printf("doCL := %d\n", doCL);
  }
 
  /* If only file name for parameter constraints was given, then write one
     with default contents, and exit */
  if(limfile[0] && !ttacfile[0]) {
    /* Check that initial value file does not exist */
    if(strcasecmp(limfile, "stdout")!=0 && access(limfile, 0) != -1) {
      fprintf(stderr, "Error: parameter constraint file %s exists.\n", limfile);
      return(9);
    }
    if(verbose>1 && strcasecmp(limfile, "stdout")!=0) 
      printf("writing parameter constraints file\n");
    /* Create parameter file */
    IFT ift; iftInit(&ift);
    iftPutDouble(&ift, "R1_lower", def_pmin[0], NULL, 0);
    iftPutDouble(&ift, "R1_upper", def_pmax[0], NULL, 0);
    iftPutDouble(&ift, "k2_lower", def_pmin[1], NULL, 0);
    iftPutDouble(&ift, "k2_upper", def_pmax[1], NULL, 0);
    iftPutDouble(&ift, "k3_lower", def_pmin[2], NULL, 0);
    iftPutDouble(&ift, "k3_upper", def_pmax[2], NULL, 0);
    ret=iftWrite(&ift, limfile, 0);
    if(ret) {
      fprintf(stderr, "Error in writing '%s': %s\n", limfile, ift.status);
      iftEmpty(&ift); return(9);
    }
    if(strcasecmp(limfile, "stdout")!=0)
      fprintf(stdout, "Parameter file %s with initial values written.\n", limfile);
    iftEmpty(&ift); return(0);
  }
 
 
  /* Did we get all the information that we need? */
  if(!resfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
 
  /*
   *  Read model parameter upper and lower limits
   *  if file for that was given
   */
  if(limfile[0]) {
    if(verbose>1) printf("reading %s\n", limfile);
    IFT ift; iftInit(&ift);
    double v;
    ret=iftRead(&ift, limfile, 1, 0);
    if(ret) {
      fprintf(stderr, "Error in reading '%s': %s\n", limfile, ift.status);
      return(9);
    }
    if(verbose>10) iftWrite(&ift, "stdout", 0);
    int n=0;
    if(iftGetDoubleValue(&ift, 0, "R1_lower", &v, 0)>=0) {def_pmin[0]=v; n++;}
    if(iftGetDoubleValue(&ift, 0, "R1_upper", &v, 0)>=0) {def_pmax[0]=v; n++;}
    if(iftGetDoubleValue(&ift, 0, "k2_lower", &v, 0)>=0) {def_pmin[1]=v; n++;}
    if(iftGetDoubleValue(&ift, 0, "k2_upper", &v, 0)>=0) {def_pmax[1]=v; n++;}
    if(iftGetDoubleValue(&ift, 0, "k3_lower", &v, 0)>=0) {def_pmin[2]=v; n++;}
    if(iftGetDoubleValue(&ift, 0, "k3_upper", &v, 0)>=0) {def_pmax[2]=v; n++;}
    iftEmpty(&ift);
    if(n==0) {fprintf(stderr, "Error: invalid parameter file.\n"); return(9);}
  }
  /* Check that these are ok */
  n=0; ret=0;
  for(int pi=0; pi<parNr; pi++) {
    if(verbose>3) printf(" %d %g %g\n", pi+1, def_pmin[pi], def_pmax[pi]);
    if(def_pmin[pi]<0.0) ret++;
    if(def_pmax[pi]<def_pmin[pi]) ret++;
    if(def_pmax[pi]>def_pmin[pi]) n++;
    if(verbose>3 && ret>0) printf("   -> invalid\n");
  }
  if(ret!=0) {
    fprintf(stderr, "Error: invalid parameter constraints.\n");
    return(9);
  }
  if(n==0) {
    fprintf(stderr, "Error: no model parameters left free for fitting.\n");
    return(9);
  }
  if(verbose>1) {
    printf("Parameter constraints:\n");
    for(int pi=0; pi<parNr; pi++) {
      printf("def_pmin[%d] := %g\n", pi+1, def_pmin[pi]);
      printf("def_pmax[%d] := %g\n", pi+1, def_pmax[pi]);
    }
  }
 
 
 
  /*
   *  Read TTAC file
   */
  if(verbose>1) printf("reading %s\n", ttacfile);
  DFT dft; dftInit(&dft);
  if(dftRead(ttacfile, &dft)) {
    fprintf(stderr, "Error in reading '%s': %s\n", ttacfile, dfterrmsg);
    return(2);
  }
  /* Check for NA's */
  if(dft_nr_of_NA(&dft) > 0) {
    fprintf(stderr, "Error: missing sample(s) in %s\n", ttacfile);
    dftEmpty(&dft); return(2);
  }
  /* Sort the data by increasing sample times */
  dftSortByFrame(&dft);
  /* Set time unit to min */
  ret=dftTimeunitConversion(&dft, TUNIT_MIN);
  if(ret) fprintf(stderr, "Warning: check that regional data times are in minutes.\n");
  /* Remove frame overlaps and gaps */
  if(dft.timetype==DFT_TIME_STARTEND) {
    if(verbose>2) fprintf(stdout, "checking frame overlap in %s\n", ttacfile);
    ret=dftDeleteFrameOverlap(&dft);
    if(ret) {
      fprintf(stderr, "Error: %s has overlapping frame times.\n", ttacfile);
      dftEmpty(&dft); return(2);
    }
  }
  /* Set fit duration */
  int first, last;
  double starttime, endtime;
  starttime=0.0; endtime=fitdur;
  fitframeNr=fittime_from_dft(&dft, &starttime, &endtime, &first, &last, verbose-2);
  if(fitframeNr<5) {
    fprintf(stderr, "Error: too few data points for a decent fit.\n");
    dftEmpty(&dft); return(2);
  }
  if(verbose>2) {
    printf("dft.frameNr := %d\n", dft.frameNr);
    printf("starttime := %g\n", starttime);
    printf("endtime := %g\n", endtime);
    printf("first := %d\n", first);
    printf("last := %d\n", last);
    printf("fitframeNr := %d\n", fitframeNr);
  }
  fitdur=endtime;
  /* Check that there is not any significant delay in the beginning of the data */
  if(dft.timetype==DFT_TIME_STARTEND) {
    if(dft.x1[0]>0.45) {
      fprintf(stderr, "Error: TACs must start at time zero.\n");
      dftEmpty(&dft); return(2);
    }
    if(dft.x1[0]>0.0833333) {
      fprintf(stderr, "Warning: TACs should start at time zero.\n");
    }
  }
  if(verbose>2) printf("Tissue calibration unit := %s\n", dft.unit);
 
  /* Add data weights, if requested */
  if(weights==1) {
    dft.isweight=0; 
    for(int i=0; i<dft.frameNr; i++) dft.w[i]=1.0;
  } else if(weights==2) {
    if(dftWeightByFreq(&dft)!=0) {
      fprintf(stderr, "Error: cannot set data weights.\n");
      dftEmpty(&dft); return(2);
    }
  } else if(dft.isweight==0) {
    fprintf(stderr, "Warning: data is not weighted.\n");
  }
  /* Print the weights */
  if(verbose>2) {
    fprintf(stdout, "common_data_weights := %g", dft.w[0]);
    for(int i=1; i<dft.frameNr; i++) fprintf(stdout, ", %g", dft.w[i]);
    fprintf(stdout, "\n");
  }
 
 
  /*
   *  Read reference TAC
   */
  if(verbose>1) printf("\nreading reference\n");
  int inputtype=-1, ref=-1;
  ret=dftReadReference(&dft, rtacfile, &inputtype, &ref, tmp, verbose-1);
  if(ret<=0) {
    fprintf(stderr, "Error in reading reference input: %s\n", tmp);
    dftEmpty(&dft); if(verbose>1) printf("ret := %d\n", ret);
    return(3);
  }
  if(ret>1) {
    fprintf(stderr, "Warning: several reference regions found: %s selected.\n", dft.voi[ref].name);
  } else if(verbose>1) printf("reference_region := %s\n", dft.voi[ref].name);
  if(verbose>2) printf("inputtype := %d\n", inputtype);
 
  /* Allocate an extra TAC for the bootstrap */
  int bsi=-1;
  if(doBootstrap) {
    ret=dftAddmem(&dft, 1); if(ret) {
      fprintf(stderr, "Error: cannot allocate more memory.\n");
      dftEmpty(&dft); return(4);
    }
    bsi=dft.voiNr;
    strcpy(dft.voi[bsi].voiname, "BS");
    strcpy(dft.voi[bsi].name, "BS");
  }
 
 
 
  /*
   *  Prepare the room for results
   */
  if(verbose>1) printf("initializing result data\n");
  RES res; resInit(&res);
  ret=res_allocate_with_dft(&res, &dft); if(ret!=0) {
    fprintf(stderr, "Error: cannot set-up memory for results.\n");
    dftEmpty(&dft); return(4);
  }
  /* Copy titles & filenames */
  tpcProgramName(argv[0], 1, 1, res.program, 256);
  strcpy(res.datafile, ttacfile);
  if(inputtype!=5 && rtacfile[0]) strcpy(res.reffile, rtacfile);
  if(ref>=0) strcpy(res.refroi, dft.voi[ref].name);
  strcpy(res.fitmethod, "TGO");
  /* Constants */
  res.isweight=dft.isweight;
  /* Set data range */
  sprintf(res.datarange, "%g - %g %s", 0.0, fitdur, petTunit(dft.timeunit));
  res.datanr=fitframeNr;
  /* Set current time to results */
  res.time=time(NULL);
  /* Set parameter number, including also the extra "parameters"
     and the parameter names and units */
  res.parNr=4;
  {
    int pi;
    pi=0; strcpy(res.parname[pi], "R1"); strcpy(res.parunit[pi], "");
    pi++; strcpy(res.parname[pi], "k2"); strcpy(res.parunit[pi], "1/min");
    pi++; strcpy(res.parname[pi], "k3"); strcpy(res.parunit[pi], "1/min");
    pi++; strcpy(res.parname[pi], "WSS"); strcpy(res.parunit[pi], "");
  }
 
 
 
  /*
   *  Fit one VOI at a time
   */
  if(verbose>0) printf("\nfitting...\n");
  int tgoNr=0, neighNr=0, iterNr=0;
  double wss;
  /* Set common data pointers */
  t=dft.x; cr=dft.voi[ref].y;
  /* Fit model to one TAC at a time */
  for(int ri=0; ri<dft.voiNr; ri++) if(ri!=ref) {
 
    if(verbose>1) printf("Region %d %s\n", ri+1, dft.voi[ri].name);
    /* Set data pointers */
    tis=dft.voi[ri].y; ct=dft.voi[ri].y2; w=dft.w;
    double *p=res.voi[ri].parameter;
 
    /* Set common parameter constraints */
    for(int pi=0; pi<parNr; pi++) {pmin[pi]=def_pmin[pi]; pmax[pi]=def_pmax[pi];}
    if(verbose>30) {
      printf("Parameter constraints:\n");
      for(int pi=0; pi<parNr; pi++) printf("  %10.3E - %10.3E\n", pmin[pi], pmax[pi]);
    }
 
    /* Fit */
    if(verbose>2) printf("  fitting curve...\n");
    TGO_LOCAL_INSIDE=0;
    TGO_SQUARED_TRANSF=0;
    tgoNr=220;
    neighNr=20;
    ret=tgo(pmin, pmax, rrtmFunc, NULL, parNr, neighNr, &wss, p, tgoNr, iterNr, verbose-8);
    if(ret>0) {
      fprintf(stderr, "Error in optimization (%d).\n", ret);
      dftEmpty(&dft); resEmpty(&res); return(6);
    }
    if(verbose>3) {
      for(int pi=0; pi<parNr; pi++) printf(" %g", p[pi]);
      printf(" -> WSS=%g\n", p[parNr]); fflush(stdout);
    }
    /* Correct fitted parameters to match constraints like inside the function */
    (void)modelCheckParameters(parNr, pmin, pmax, p, p, NULL);
    p[parNr]=wss=wss_wo_penalty;
    if(verbose>2) printf("wss := %g\nfitframeNr := %d\n", wss, fitframeNr);
 
    /* Bootstrap */
    if(doBootstrap) {
      if(verbose>2) printf("  bootstrapping...\n");
      /* bootstrap changes measured and simulated data, therefore use copies */
      tis=dft.voi[bsi].y; ct=dft.voi[bsi].y2;
      if(doSD) sd=res.voi[ri].sd; else sd=NULL;
      if(doCL) {cl1=res.voi[ri].cl1; cl2=res.voi[ri].cl2;} else cl1=cl2=NULL;
      ret=bootstrap(
        0, cl1, cl2, sd, p, pmin, pmax, fitframeNr,
        // measured and fitted original TAC, not modified
        dft.voi[ri].y, dft.voi[ri].y2,
        // tissue TAC noisy data is written to be used by objf
        tis, 
        parNr, w, rrtmFunc, tmp, verbose-5
      );
      if(ret) {
        fprintf(stderr, "Error in bootstrap: %s\n", tmp);
        for(int pi=0; pi<parNr; pi++) {
          if(doSD) sd[pi]=nan(""); 
          if(doCL) cl1[pi]=cl2[pi]=nan("");
        }
      }
      // back to what pointers were
      tis=dft.voi[ri].y; tis=dft.voi[ri].y2;
    }
 
  } /* Next VOI */
  if(verbose>0) {fprintf(stdout, "\n"); fflush(stdout);}
 
  /* Delete the reference region from the results */
  resDelete(&res, ref);
  /* If reference data file was given and it contained several TACs then delete those too */
  if(inputtype!=5) {
    for(int i=dft.voiNr-1; i>=0; i--) if(dft.voi[i].sw!=0) resDelete(&res, i);
  }
 
  /*
   *  Print results on screen
   */
  if(verbose>0) {resPrint(&res); fprintf(stdout, "\n");}
 
 
  /*
   *  Save results
   */
  if(verbose>1) printf("saving results in %s\n", resfile);
  ret=resWrite(&res, resfile, verbose-5);
  if(ret) {
    fprintf(stderr, "Error in writing '%s': %s\n", resfile, reserrmsg);
    resEmpty(&res); dftEmpty(&dft);
    return(11);
  }
  if(verbose>0) fprintf(stdout, "Model parameters written in %s\n", resfile);
 
 
 
  /*
   *  Saving and/or plotting of fitted TACs
   */
  if(svgfile[0] || fitfile[0]) {
 
    /* Create a DFT containing fitted TACs */
    char tmp[64];
    DFT dft2;
    dftInit(&dft2); ret=dftdup(&dft, &dft2);
    if(ret) {
      fprintf(stderr, "Error: cannot save fitted curves.\n");
      dftEmpty(&dft); resEmpty(&res);
      return(21);
    }
    for(int ri=0; ri<dft.voiNr; ri++) 
      if(ri!=ref) 
        for(int fi=0; fi<fitframeNr; fi++)
          dft2.voi[ri].y[fi]=dft2.voi[ri].y2[fi];
    dft2.frameNr=fitframeNr;
 
    /* Save SVG plot of fitted and original data */
    if(svgfile[0]) {
      if(verbose>1) printf("saving SVG plot\n");
      sprintf(tmp, "RRTM fit ");
      if(strlen(dft.studynr)>0) strcat(tmp, dft.studynr);
      ret=plot_fitrange_svg(&dft, &dft2, tmp, 0.0, 1.02*dft.x[fitframeNr-1],
                            0.0, nan(""), svgfile, verbose-8);
      if(ret) {
        fprintf(stderr, "Error (%d) in writing '%s'.\n", ret, svgfile);
        dftEmpty(&dft2); dftEmpty(&dft); resEmpty(&res);
        return(30+ret);
      }
      if(verbose>0) printf("Plots written in %s\n", svgfile);
    }
 
    /* Delete reference region(s) from the data, unless it already existed in data */
    if(inputtype!=5) {
      for(int i=dft2.voiNr-1; i>=0; i--) if(dft2.voi[i].sw!=0) dftDelete(&dft2, i);
    }
 
    /* Save fitted TACs */
    if(fitfile[0]) {
      if(verbose>1) printf("saving fitted curves\n");
      tpcProgramName(argv[0], 1, 0, tmp, 128);
      sprintf(dft2.comments, "# program := %s\n", tmp);
      if(dftWrite(&dft2, fitfile)) {
        fprintf(stderr, "Error in writing '%s': %s\n", fitfile, dfterrmsg);
        dftEmpty(&dft2); dftEmpty(&dft); resEmpty(&res);
        return(22);
      }
      if(verbose>0) printf("Fitted TACs written in %s\n", fitfile);
    }
 
    dftEmpty(&dft2);
  }
 
 
  resEmpty(&res);
  dftEmpty(&dft);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************
 *
 *  Functions to be minimized
 *
 *****************************************************************************/
double rrtmFunc(int parNr, double *p, void *fdata)
{
  int ret;
  double R1, k2, k3, d, wss=0.0;
  double pa[MAX_PARAMETERS], penalty=1.0;
 
 
  /* Check parameters against the constraints */
  ret=modelCheckParameters(parNr, pmin, pmax, p, pa, &penalty);
  if(fdata) {}
  /* Get parameters */
  R1=pa[0]; k2=pa[1]; k3=pa[2];
 
  /* Simulate the tissue PET TAC */
  ret=simRTCM(t, cr, fitframeNr, R1, k2, k3, 0.0, ct, NULL, NULL);
  if(ret) {
    fprintf(stderr, "  error %d in simulation\n", ret);
    return(nan(""));
  }
 
  /* Calculate error */
  for(int i=0; i<fitframeNr; i++) if(w[i]>0.0) {
    d=ct[i]-tis[i]; 
    wss+=w[i]*d*d;
  }
  wss_wo_penalty=wss;
  wss*=penalty;
  if(0) printf("R1=%g  k2=%g  k3=%g  => %g\n", R1, k2, k3, wss);
 
  return(wss);
}
/*****************************************************************************/
 
/*****************************************************************************/