fit_trtm.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 trtmFunc(int parNr, double *p, void*);
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "NLLSQ estimation of R1 (=K1/K1'), k2', and k3 applying the transport-limited",
  "reference tissue compartment model, TRTM (1). This model is based on",
  "the reference tissue compartment model (2), but here it is assumed that",
  "metabolism of tracer is irreversible (k4=0) during the PET scanning, and that",
  "in (positive) reference tissue k3'>>k2', and thus the uptake in it is limited",
  "only by transport into tissue (3,4).",
  " ",
  "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:",
  " -rk2=<<value>|mean|median>",
  "     Constrain k2' to specified <value>, or to mean or median of regional",
  "     k2' values excluding reference region(s)",
  " -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.tac contains regions-of-interest and positive reference",
  "region, with name 'putam'. The whole time range is used in the fit.",
  "Fitted TACs are plotted in SVG format in file a789trtm.svg.",
  "     @P -svg=a789trtm.svg a789.tac putam 999 a789.res",
  " ",
  "Example 2: Reference region TAC is in a separate file, a789ref.tac.",
  "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.tac a789ref.tac 60 a789.res",
  " ",
  "Example 3a: Create a file containing default parameter limits:",
  "     @P -lim=trtm.lim",
  " ",
  "Example 3b: Apply user-defined parameter constraints specified in trtm.lim:",
  "     @P -lim=trtm.lim a789.tac 'putam mean' 999 a789.tac",
  " ",
  "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. Herholz K, Lercher M, Wienhard K, Bauer B, Lenz O, Heiss W-D.",
  "   PET measurement of cerebral acetylcholine esterase activity without",
  "   blood sampling. Eur J Nucl Med 2001;28:472-477.",
  "4. Nagatsuka S, Fukushi K, Shinotoh H, Namba H, Iyo M, Tanaka N, Aotsuka A,",
  "   Ota T, Tanada S, Irie T. Kinetic analysis of [11C]MP4A using a high-",
  "   radioactivity brain region that represents an integrated input function",
  "   for measurement of cerebral acetylcholinesterase activity without",
  "   arterial blood sampling. J Cereb Blood Flow Metab 2001; 21: 1354-1366.",
  " ",
  "See also: fitk3, tacweigh, fit_rrtm, lhtrtm, sim_rtcm, rescoll",
  " ",
  "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];
  int      fix_rk2=0;
  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]=10.0;
  /* rk2 */ def_pmin[1]=0.000001;  def_pmax[1]=1.0;
  /* k3  */ def_pmin[2]=0.0;       def_pmax[2]=10.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;
    } else if(strncasecmp(cptr, "RK2=", 4)==0 && strlen(cptr)>4) {
      cptr+=4;
      if(strcasecmp(cptr, "MEAN")==0) {fix_rk2=2; continue;}
      else if(strncasecmp(cptr, "AVERAGE", 2)==0) {fix_rk2=2; continue;}
      if(strcasecmp(cptr, "MEDIAN")==0) {fix_rk2=3; continue;}
      fix_rk2=1; def_pmin[1]=def_pmax[1]=atof_dpi(cptr);
      if(def_pmin[1]<=0.0) {
        fprintf(stderr, "Error: invalid value with option '%s'.\n", argv[ai]);
        return(1);
      }
      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);
    printf("fix_rk2 := %d\n", fix_rk2);
    if(fix_rk2==1) printf("fixed_rk2 := %g\n", def_pmin[1]);
  }
 
  /* 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, "rk2_lower", def_pmin[1], NULL, 0);
    iftPutDouble(&ift, "rk2_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, "rk2_lower", &v, 0)>=0) {n++; if(fix_rk2!=1) def_pmin[1]=v;}
    if(iftGetDoubleValue(&ift, 0, "rk2_upper", &v, 0)>=0) {n++; if(fix_rk2!=1) def_pmax[1]=v;}
    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], "rk2"); 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], "");
  }
 
 
  /*
   *  Determine average and median k2' by fitting, if regional k2' will be fixed to it
   */
  int tgoNr=0, neighNr=0, iterNr=0;
  double wss;
  if(fix_rk2 >= 2) {
    double rk2median, rk2mean, rk2sd, rk2list[res.voiNr], v;
 
    if(verbose>0) {
      if(fix_rk2==2) printf("\nfitting k2' for fixing to its regional mean...\n");
      else printf("\nfitting k2' for fixing to its regional median...\n");
    }
 
    /* Set common data pointers */
    t=dft.x;
    cr=dft.voi[ref].y;
 
    /* Fit regions */
    for(int ri=0; ri<dft.voiNr; ri++) if(dft.voi[ri].sw==0) {
      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 parameter constraints */
      for(int pi=0; pi<parNr; pi++) {pmin[pi]=def_pmin[pi];  pmax[pi]=def_pmax[pi];}
      if(verbose>30) {
        printf("  constraints :=");
        for(int pi=0; pi<parNr; pi++) printf(" [%g,%g]", pmin[pi], pmax[pi]);
        printf("\n");
      }
      /* Fit */
      if(verbose>2) printf("  fitting curve...\n");
      TGO_LOCAL_INSIDE=0;
      TGO_SQUARED_TRANSF=0;
      tgoNr=220;
      neighNr=20;
      ret=tgo(pmin, pmax, trtmFunc, 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);
    } // next region
 
    /* Calculate average and median k2' */
    for(int ri=0, n=0; ri<res.voiNr; ri++) if(dft.voi[ri].sw==0) {
      /* check that k2 is reasonable */
      if(isnan(res.voi[ri].parameter[1]) || res.voi[ri].parameter[1]<=0.010) continue;
      rk2list[n++]=res.voi[ri].parameter[1];
    }
    rk2mean=dmean(rk2list, n, &rk2sd);
    rk2median=dmedian(rk2list, n);
    if(verbose>1) {
      if(verbose>2) printf("n_rk2 := %d\n", n);
      printf("mean_rk2 := %6.4f +- %6.4f\n", rk2mean, rk2sd);
      printf("median_rk2 := %6.4f\n", rk2median);
    }
 
    /* Fix rk2' */
    if(fix_rk2==2) v=rk2mean; else v=rk2median;
    if(v<def_pmin[1] || v>def_pmax[1]) {
      fprintf(stderr, "Error: k2' could not be fixed to regional mean or median\n");
      dftEmpty(&dft); resEmpty(&res); return(7);
    }
    def_pmin[1]=def_pmax[1]=v;
    if(verbose>1) printf("fixed_rk2 := %g\n", def_pmin[1]);
  }
 
 
  /*
   *  Fit one VOI at a time
   */
  if(verbose>0) printf("\nfitting...\n");
  /* 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;
    if(fix_rk2==1) {
      tgoNr=220;
      neighNr=20;
    } else {
      tgoNr=180;
      neighNr=16;
    }
    ret=tgo(pmin, pmax, trtmFunc, 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, trtmFunc, 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, "TRTM 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 trtmFunc(int parNr, double *p, void *fdata)
{
  int ret;
  double R1, k2, rk2, 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]; rk2=pa[1]; k3=pa[2];
  k2=R1*rk2;
 
  /* Simulate the tissue PET TAC */
  ret=simTRTM(t, cr, fitframeNr, R1, k2, k3, ct);
  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);
}
/*****************************************************************************/
 
/*****************************************************************************/