fit_h2o.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"
/*****************************************************************************/
 
/*****************************************************************************/
int parNr=4;
DFT btac, ttac;
double *petmeas, *petsim;
double fVa=-1.0, fK1k2=-1.0, fDelay;
double pmin[MAX_PARAMETERS], pmax[MAX_PARAMETERS];
int fitframeNr;
double wss_wo_penalty=0.0;
/* pointer to data, not allocated */
double *ca, *cs, *cs2, *simdc, *petdc, *weight, *cpet, *csim, *wres;
/* Local functions */
double waterFunc(int parNr, double *p, void*);
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Non-linear fitting of one-tissue compartmental water model to regional",
  "dynamic PET [O-15]H2O study data.",
  "The model parameters are blood flow (F), partition coefficient",
  "(pWater, K1/k2), arterial blood volume (Va) and time delay (delayT) between",
  "tissue and blood time-activity curves (TTAC and BTAC, respectively).",
  "Venous radioactivity is assumed to be the same as the tissue concentration.",
  "Extraction factor of water is assumed to be 1.0.",
  " ",
  "    ______        ______    ",
  "   |      |  K1  |      |   ",
  "   |  Ca  | ---> |  C1  |   ",
  "   |______|      |______|   ",
  "                    | k2    ",
  "                    V       ",
  " ",
  "The blood flow obtained using PET [O-15]H2O techniques reflects tissue",
  "perfusion, since the method is dependent on the tissue exchange of",
  "labelled water. Non-nutritive flow (blood flowing through arteriovenous",
  "shunts) is not measured (Lammertsma & Jones, 1983).",
  " ",
  "Usage: @P [Options] btacfile ttacfile endtime resultfile",
  " ",
  "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.",
  "     Without filename the default values are printed on screen.",
  "     Parameter can be fixed to a certain value by setting its",
  "     lower and upper limit to that value.",
  " -ml or -dl",
  "     In results the units of F and Va will be given per mL or per dL,",
  "     respectively. By default, units are per dL.",
  " -fpt",
  "     Blood flow (perfusion) is reported per tissue volume (PET volume minus",
  "     vascular volume), that is, F=K1. By default, perfusion is reported",
  "     per PET volume, that is, F=(1-Va)*K1.",
  " -SD[=<y|N>]",
  "     Standard deviations are calculated and saved in results (Y, default),",
  "     or not calculated (n).",
  "     Program runs a lot faster if SD and CL are not calculated.",
  " -CL[=<y|N>]",
  "     95% Confidence limits are calculated and saved in results (y), or",
  "     not calculated (N, default).",
  " -Va=<Va(%)>",
  "     Enter a fixed Va; fitted by default.",
  " -Delay=<Time delay (s)>",
  "     Enter a fixed time delay; fitted by default.",
  "     Positive time delay means that BTAC is moved forward in time.",
  " -pH2O=<Partition coefficient for water>",
  "     Enter a fixed pH2O; fitted by default.",
  " -fit=<Filename>",
  "     Fitted regional TACs are written in DFT format.",
  " -input=<Filename>",
  "     Input TAC sample times are corrected by the median of fitted time",
  "     delay values and saved; resulting input file can be used with imgflow,",
  "     or as input to this program to have common time delay for all regions.",
  " -resid=<Filename>",
  "     Weighted residual curves are written in DFT format.",
  " -svg=<Filename>",
  "     Fitted and measured TACs are plotted in specified SVG file.",
  " -k2",
  "     Save k2 in result file.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example:",
  "     @P uo1234bl.kbq uo1234.dft 999 uo1234f.res",
  " ",
  "See also: b2t_h2o, imgflow, arlkup, water_input, fitk2, bfmh2o",
  " ",
  "Keywords: TAC, modelling, perfusion, blood flow, 1TCM",
  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;
  double     def_pmin[MAX_PARAMETERS], def_pmax[MAX_PARAMETERS];
  char       btacfile[FILENAME_MAX], ttacfile[FILENAME_MAX], 
             limfile[FILENAME_MAX], resfile[FILENAME_MAX], 
             fitfile[FILENAME_MAX], residfile[FILENAME_MAX],
             inputfile[FILENAME_MAX], svgfile[FILENAME_MAX];
  int        doBootstrap=0, doSD=0, doCL=0; // 0=no, 1=yes
  int        per_dl=1; // 0 or 1
  int        flow_per_tissue=0; // 0: f=(1-Va)K1 ; 1: f=K1
  int        save_k2=0;
  double     fitdur=-1.0;
  int        fittedparNr=0;
 
  IFT        ift;
  RES        res;
  char      *cptr, tmp[FILENAME_MAX];
  int        times_changed=0;
  int        sim_index, simpet_index, bs_index;
  double     lambda=0.0;
  double     wss;
 
 
  /* Set lambda for O-15 */
  lambda=hl2lambda(HL_O15*60.0);
 
  /* Set parameter initial values and constraints */
  fDelay=nan("");
  /* Flow  */ def_pmin[0]=0.0;       def_pmax[0]=600.0;
  /* pH2O  */ def_pmin[1]=0.30;      def_pmax[1]=1.0;
  /* Va    */ def_pmin[2]=0.0;       def_pmax[2]=25.0;
  /* Delay */ def_pmin[3]=-60.0;     def_pmax[3]=+60.0;
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  limfile[0]=btacfile[0]=ttacfile[0]=resfile[0]=inputfile[0]=(char)0;
  svgfile[0]=fitfile[0]=residfile[0]=(char)0;
  iftInit(&ift);
  dftInit(&btac); dftInit(&ttac);
  resInit(&res);
  /* 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(strncasecmp(cptr, "I=", 2)==0 && strlen(cptr)>2) { // deprecated
      strlcpy(limfile, cptr+2, FILENAME_MAX); continue;
    } else if(strcasecmp(cptr, "LIM")==0) {
      strcpy(limfile, "stdout"); continue;
    } else if(strcasecmp(cptr, "I")==0) { // deprecated
      strcpy(limfile, "stdout"); continue;
    } else if(strcasecmp(cptr, "DL")==0) {
      per_dl=1; continue;
    } else if(strcasecmp(cptr, "ML")==0) {
      per_dl=0; continue;
    } else if(strcasecmp(cptr, "FPT")==0) {
      flow_per_tissue=1; continue;
    } else if(strcasecmp(cptr, "K2")==0) {
      save_k2=1; continue;
    } if(strncasecmp(cptr, "INPUT=", 6)==0) {
      strlcpy(inputfile, cptr+6, FILENAME_MAX); 
      if(strlen(inputfile)>0) continue;
    } else if(strncasecmp(cptr, "SVG=", 4)==0) {
      strlcpy(svgfile, cptr+4, FILENAME_MAX); 
      if(strlen(svgfile)>0) continue;
    } else if(strncasecmp(cptr, "FIT=", 4)==0) {
      strlcpy(fitfile, cptr+4, FILENAME_MAX); 
      if(strlen(fitfile)>0) continue;
    } else if(strncasecmp(cptr, "RESID=", 6)==0) {
      strlcpy(residfile, cptr+6, FILENAME_MAX); 
      if(strlen(residfile)>0) continue;
    } else if((strncasecmp(cptr, "Va=", 3)==0 || strncasecmp(cptr, "Vb=", 3)==0) && strlen(cptr)>3) {
      fVa=atof_dpi(cptr+3);
      if(fVa>=0.0 && fVa<100.0) {
        def_pmin[2]=def_pmax[2]=fVa;
        if(fVa>0.0 && fVa<0.5) 
          fprintf(stderr, "Warning: Va was set to %g%%\n", fVa);
        continue;
      }
    } else if(strncasecmp(cptr, "pH2O=", 5)==0 && strlen(cptr)>5) {
      fK1k2=atof_dpi(cptr+5); if(fK1k2>0.0) continue;
    } else if(strncasecmp(cptr, "delay=", 6)==0 && strlen(cptr)>6) {
      fDelay=atof_dpi(cptr+6);
      if(fDelay>-120.0 && fDelay<120.0) {
        if(fDelay!=0.0 && fabs(fDelay)<1.0)
          fprintf(stderr, "Warning: Delay was set to %g sec\n", fDelay);
        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(btacfile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {strlcpy(ttacfile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {
    if(atof_with_check(argv[ai], &fitdur) || fitdur<0.0) {
      fprintf(stderr, "Error: invalid fit time '%s'.\n", argv[ai]); return(1);}
    ai++;
  }
  if(ai<argc) {strlcpy(resfile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {fprintf(stderr, "Error: too many arguments: '%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("limfile := %s\n", limfile);
    printf("btacfile := %s\n", btacfile);
    printf("ttacfile := %s\n", ttacfile);
    printf("resfile := %s\n", resfile);
    printf("fitfile := %s\n", fitfile);
    printf("svgfile := %s\n", svgfile);
    printf("residfile := %s\n", residfile);
    printf("inputfile := %s\n", inputfile);
    printf("per_dl := %d\n", per_dl);
    printf("flow_per_tissue := %d\n", flow_per_tissue);
    printf("doBootstrap := %d\n", doBootstrap);
    printf("doSD := %d\n", doSD);
    printf("doCL := %d\n", doCL);
    if(fitdur>=0.0) printf("required_fitdur := %g\n", fitdur);
    if(fVa>=0.0) printf("required_fVa := %g\n", fVa);
    if(fK1k2>=0.0) printf("required_fK1k2 := %g\n", fK1k2);
    if(!isnan(fDelay)) printf("required_fDelay := %g\n", fDelay);
    printf("save_k2 := %d\n", save_k2);
    printf("lambda := %g [1/s]\n", lambda);
    fflush(stdout);
  }
 
  /* If only filename for initial values was given, then write one
     with default contents, and exit */
  if(limfile[0] && !btacfile[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) printf("writing parameter constraints file\n");
    /* Create parameter file */
    iftPutDouble(&ift, "K1_lower", def_pmin[0], NULL, 0);
    iftPutDouble(&ift, "K1_upper", def_pmax[0], NULL, 0);
    iftPutDouble(&ift, "K1k2_lower", def_pmin[1], NULL, 0);
    iftPutDouble(&ift, "K1k2_upper", def_pmax[1], NULL, 0);
    iftPutDouble(&ift, "Va_lower", def_pmin[2], NULL, 0);
    iftPutDouble(&ift, "Va_upper", def_pmax[2], NULL, 0);
    iftPutDouble(&ift, "Delay_lower", def_pmin[3], NULL, 0);
    iftPutDouble(&ift, "Delay_upper", def_pmax[3], NULL, 0);
    if(iftWrite(&ift, limfile, 0)) {
      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(fitdur==0) fitdur=1.0E+100;
  else if(fitdur<0) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  if(!resfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
 
  /*
   *  Read model parameter initial values and upper and lower limits
   *  if file for that was given
   */
  if(limfile[0]) {
    double v;
    if(verbose>1) printf("reading %s\n", limfile);
    if(iftRead(&ift, limfile, 1, 0)) {
      fprintf(stderr, "Error in reading '%s': %s\n", limfile, ift.status);
      return(9);
    }
    if(verbose>10) iftWrite(&ift, "stdout", 0);
    int n=0;
    /* K1 */
    if(iftGetDoubleValue(&ift, 0, "K1_lower", &v, 0)>=0) {def_pmin[0]=v; n++;}
    if(iftGetDoubleValue(&ift, 0, "K1_upper", &v, 0)>=0) {def_pmax[0]=v; n++;}
    /* K1/k2 */
    if(iftGetDoubleValue(&ift, 0, "K1k2_lower", &v, 0)>=0) {def_pmin[1]=v; n++;}
    if(iftGetDoubleValue(&ift, 0, "K1k2_upper", &v, 0)>=0) {def_pmax[1]=v; n++;}
    /* Va */
    if(iftGetDoubleValue(&ift, 0, "Va_lower", &v, 0)>=0) {def_pmin[2]=v; n++;}
    if(iftGetDoubleValue(&ift, 0, "Va_upper", &v, 0)>=0) {def_pmax[2]=v; n++;}
    /* Delay */
    if(iftGetDoubleValue(&ift, 0, "Delay_lower", &v, 0)>=0) {def_pmin[3]=v; n++;}
    if(iftGetDoubleValue(&ift, 0, "Delay_upper", &v, 0)>=0) {def_pmax[3]=v; n++;}
    iftEmpty(&ift);
    if(n==0) {fprintf(stderr, "Error: invalid parameter file.\n"); return(9);}
  }
  /* Check that these limits are ok */
  {
    int n=0, pi, ret=0;
    for(pi=0; pi<parNr; pi++) {
      if(def_pmin[pi]<0.0 && pi!=3) ret++;
      if(def_pmax[pi]<def_pmin[pi]) ret++;
      if(def_pmax[pi]>def_pmin[pi]) n++;
    }
    if(ret) {
      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);
    }
  }
  /* Fixed/fitted Va */
  if(fVa>=0.0) def_pmin[2]=def_pmax[2]=fVa;
  else if(def_pmin[2]==def_pmax[2]) fVa=def_pmin[2];
  if(verbose>1) {
    if(fVa>=0.0) printf("fVa := %g\n", fVa);
  }
  /* Fixed/fitted pH2O */
  if(fK1k2>=0.0) def_pmin[1]=def_pmax[1]=fK1k2;
  else if(def_pmin[1]==def_pmax[1]) fK1k2=def_pmin[1];
  if(verbose>1) {
    if(fK1k2>=0.0) printf("fK1k2 := %g\n", fK1k2);
  }
  /* Fixed/fitted time delay */
  if(!isnan(fDelay)) def_pmin[3]=def_pmax[3]=fDelay;
  else if(def_pmin[3]==def_pmax[3]) fDelay=def_pmin[3];
  if(verbose>1) {
    if(!isnan(fDelay)) printf("fDelay := %g\n", fDelay);
  }
  /* Convert parameter settings */
  /* F mL/min/dL -> mL/sec/mL */ def_pmin[0]/=6000.; def_pmax[0]/=6000.;
  /* Va mL/dL    -> mL/mL */     def_pmin[2]/=100.;  def_pmax[2]/=100.;
  if(verbose>1) {
    fflush(stdout); 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]);
    }
    fflush(stdout);
  }
 
  /*
   *  Read tissue and input data
   */
  if(verbose>1) printf("reading tissue and input data\n");
  fitdur/=60.0; // sec -> min
  /* Read blofile twice, the copy is used to store simulated TAC */
  if(dftReadModelingData(ttacfile, btacfile, btacfile, NULL, &fitdur, 
                          &fitframeNr, &ttac, &btac, stdout, verbose-2, tmp)!=0)
  {
    fprintf(stderr, "Error: %s\n", tmp);
    return(2);
  }
  if(fitframeNr<4 || btac.frameNr<4) {
    if(ttac.timeunit==TUNIT_UNKNOWN || btac.timeunit==TUNIT_UNKNOWN)
      fprintf(stdout, "Error: check the time units in data files.\n");
    else
      fprintf(stderr, "Error: too few samples in specified fit duration.\n");
    dftEmpty(&btac); dftEmpty(&ttac); return(3);
  }
  sim_index=btac.voiNr-1;
  /* Prepare place for simulated tissue TAC and for bootstrap;
     do this before changing frameNr */
  if(dftAddmem(&ttac, 2)) {
    fprintf(stderr, "Error: cannot setup memory for simulation.\n");
    dftEmpty(&btac); dftEmpty(&ttac); return(3);
  }
  simpet_index=ttac.voiNr; bs_index=ttac.voiNr+1; 
  strcpy(ttac.voi[simpet_index].voiname, "Sim"); 
  strcpy(ttac.voi[simpet_index].name, "Sim");
  strcpy(ttac.voi[bs_index].voiname, "BS");
  strcpy(ttac.voi[bs_index].name, "BS");
  /* Times should be in minutes, change them to sec */
  if(ttac.timeunit==TUNIT_MIN) {
    if(verbose>1) printf("converting data times to sec\n");
    dftMin2sec(&ttac); times_changed=1; dftMin2sec(&btac);
    fitdur*=60.; // min -> sec
  } else {
    if(verbose>0) {
      printf("BTAC timeunit := %s\n", petTunit(btac.timeunit));
      printf("TTAC timeunit := %s\n", petTunit(ttac.timeunit));
    }
    if(verbose>1) printf("assuming that data times are in min.\n");
    ttac.timeunit=btac.timeunit=TUNIT_MIN;
    dftMin2sec(&ttac); times_changed=1; dftMin2sec(&btac);
    fitdur*=60.; // min -> sec
  }
  /* Check that there is not any significant delay in the beginning of data */
  if(ttac.timetype==DFT_TIME_STARTEND) {
    if(ttac.x1[0]>27.0 && ttac.voi[0].y[0]>0.0) {
      fprintf(stderr, "Error: TACs must start at time zero.\n");
      dftEmpty(&btac); dftEmpty(&ttac); return(4);
    }
    if(ttac.x1[0]>5.0 && ttac.voi[0].y[0]>0.0) {
      fprintf(stderr, "Warning: TACs should start at time zero.\n");
    }
  }
  if(btac.x[0]>90.0 && ttac.voi[0].y[0]>0.0) {
    fprintf(stderr, "Error: input TAC must start at time zero.\n");
    dftEmpty(&btac); dftEmpty(&ttac); return(4);
  }
  if(btac.x[0]>20.0 && ttac.voi[0].y[0]>0.0) {
    fprintf(stderr, "Warning: input TAC should start at time zero.\n");
  }
  /* Check weighting */
  if(ttac.isweight==0 && verbose>0)
    fprintf(stderr, "Note: data is not weighted.\n");
  /* Print the weights */
  if(verbose>1) {
    int fi=0;
    fprintf(stdout, "common_data_weights := %.3f", ttac.w[fi++]);
    for(; fi<ttac.frameNr; fi++) fprintf(stdout, ", %.3f", ttac.w[fi]);
    fprintf(stdout, "\n");
  }
 
 
 
  /*
   *  Prepare the room for results
   */
  if(verbose>1) printf("initializing result data\n");
  if(res_allocate_with_dft(&res, &ttac)!=0) {
    fprintf(stderr, "Error: cannot setup memory for results.\n");
    dftEmpty(&btac); dftEmpty(&ttac); return(6);
  }
  /* Copy titles & filenames */
  tpcProgramName(argv[0], 1, 1, res.program, 256);
  strcpy(res.bloodfile, btacfile);
  strcpy(res.datafile, ttacfile);
  strcpy(res.fitmethod, "TGO");
  /* Constants */
  if(fVa>=0.0) res.Vb=fVa;
  res.isweight=ttac.isweight;
  /* Set data range */
  sprintf(res.datarange, "%g - %g %s", 0.0, fitdur, petTunit(ttac.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=parNr+1; // WSS
  if(save_k2) res.parNr++;
  {
    int pi;
    pi=0; strcpy(res.parname[pi], "Flow"); 
    if(per_dl==0) strcpy(res.parunit[pi], "ml/(min*ml)");
    else strcpy(res.parunit[pi], "ml/(min*dl)");
    pi++; strcpy(res.parname[pi], "pWater"); strcpy(res.parunit[pi], "");
    pi++; strcpy(res.parname[pi], "Va");
    if(per_dl==0) strcpy(res.parunit[pi], "ml/ml");
    else strcpy(res.parunit[pi], "%");
    pi++; strcpy(res.parname[pi], "delayT"); strcpy(res.parunit[pi], "sec");
    if(save_k2) {pi++; strcpy(res.parname[pi], "k2"); strcpy(res.parunit[pi], "1/min");}
    pi++; strcpy(res.parname[pi], "WSS"); strcpy(res.parunit[pi], "");
  }
 
 
  /*
   *  Fit regional TACs
   */
  if(verbose>1) printf("starting regional fitting\n");
  int tgoNr=0, neighNr=5, iterNr=0;
  double *sd, *cl1, *cl2;
  if(verbose>4) {
    if(per_dl)
      fprintf(stdout, "%-20s  %-6s %-6s %-6s %-6s %-9s\n",
        "Region", "Flow", "pWater", "Va(%)", "delayT", "WSS");
    else
      fprintf(stdout, "%-20s  %-6s %-6s %-6s %-6s %-9s\n",
        "Region", "Flow", "pWater", "Va", "delayT", "WSS");
  }
  /* Calculate decay correction factors for removing decay correction from
     simulated tissue TAC (at input sample times) */
  simdc=btac.voi[sim_index].y3;
  for(int fi=0; fi<btac.frameNr; fi++)
    simdc[fi]=hlLambda2factor(-lambda, btac.x[fi], -1.0);
  /* Calculate decay correction factors for making decay correction for
     simulated tissue TAC (at PET sample times) */
  petdc=ttac.voi[simpet_index].y;
  for(int fi=0; fi<ttac.frameNr; fi++) {
    if(ttac.timetype==DFT_TIME_STARTEND)
      petdc[fi]=hlLambda2factor(lambda, ttac.x1[fi], ttac.x2[fi]-ttac.x1[fi]);
    else
      petdc[fi]=hlLambda2factor(lambda, ttac.x[fi], -1.0);
  }
  /* One region at a time */
  if(verbose>0) {
    fprintf(stdout, "fitting regional TACs: ");
    if(verbose>1) fprintf(stdout, "\n");
    fflush(stdout);
  }
  for(int ri=0; ri<ttac.voiNr; ri++) {
    int pi, ret;
    if(verbose>2) fprintf(stdout, "\n%s [%d] :\n", ttac.voi[ri].name, ri+1);
 
    /* Set data pointers */
    ca=btac.voi[0].y;
    cs=btac.voi[sim_index].y; cs2=btac.voi[sim_index].y2;
    weight=ttac.w;
    cpet=ttac.voi[ri].y; csim=ttac.voi[ri].y2; wres=ttac.voi[ri].y3;
 
    /* Set parameter constraints */
    for(pi=0; pi<parNr; pi++) {
      pmin[pi]=def_pmin[pi];  pmax[pi]=def_pmax[pi];
    }
    for(pi=fittedparNr=0; pi<parNr; pi++) if(pmax[pi]>pmin[pi]) fittedparNr++;
    if(verbose>6) {
      printf("  constraints :=");
      for(pi=0; pi<parNr; pi++) printf(" [%g,%g]", pmin[pi], pmax[pi]);
      printf("\n");
      printf("fittedparNr := %d\n", fittedparNr);
    }
 
    /* Fit */
    if(verbose>2) printf("  fitting\n");
    TGO_LOCAL_INSIDE = 1;
    TGO_SQUARED_TRANSF = 0;
    //tgoNr=0; neighNr=5;
    tgoNr=50+25*fittedparNr; /* 0, 100, 50+25 */
    neighNr=6*fittedparNr; /* 6, 10 */
    iterNr=0;
    ret=tgo(
      pmin, pmax, waterFunc, NULL, parNr, neighNr,
      &wss, res.voi[ri].parameter, tgoNr, iterNr, 
      verbose-8);
    if(ret>0) {
      fprintf(stderr, "\nError in optimization (%d).\n", ret);
      dftEmpty(&btac); dftEmpty(&ttac); resEmpty(&res); return(8);
    }
    /* Correct fitted parameters to match constraints like inside function */
    (void)modelCheckParameters(parNr, pmin, pmax, res.voi[ri].parameter,
                               res.voi[ri].parameter, NULL);
    wss=wss_wo_penalty;
    if(verbose>3) {
      for(pi=0; pi<MAX_PARAMS; pi++) printf(" %g", res.voi[ri].parameter[pi]);
      printf(" -> WSS=%g\n", wss); fflush(stdout);
    }
    res.voi[ri].parameter[res.parNr-1]=wss;
 
    /* Bootstrap */
    if(doBootstrap) {
      if(verbose>2) printf("\n  bootstrapping\n");
      /* bootstrap changes measured and simulated data, therefore use copies */
      cpet=ttac.voi[bs_index].y;
      csim=ttac.voi[bs_index].y2;
      wres=ttac.voi[bs_index].y3;
      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, res.voi[ri].parameter, pmin, pmax, fitframeNr,
        // measured original TAC, not modified
        ttac.voi[ri].y,
        // fitted TAC, not modified
        ttac.voi[ri].y2,
        // tissue TAC noisy data is written to be used by objf
        cpet,
        parNr, ttac.w, waterFunc, tmp, verbose-5
      );
      if(ret) {
        fprintf(stderr, "\nError in bootstrap: %s\n", tmp);
        for(pi=0; pi<parNr; pi++) {
          if(doSD) sd[pi]=nan(""); 
          if(doCL) cl1[pi]=cl2[pi]=nan("");
        }
      }
      // return pointers back to what they were, in case somebody adds some
      // code later after this
      cpet=ttac.voi[ri].y; csim=ttac.voi[ri].y2;
    }
 
    /* Set fitted parameters to understandable units */
    {
      double F, pWater, Va, delayT, f=60.; 
      /* k2, if requested */
      if(save_k2)
        res.voi[ri].parameter[4]=60.0*res.voi[ri].parameter[0]/res.voi[ri].parameter[1];
      /* flow */
      if(flow_per_tissue==0) f*=(1.0-res.voi[ri].parameter[2]);
      res.voi[ri].parameter[0]*=f;
      if(!isnan(res.voi[ri].cl1[0])) res.voi[ri].cl1[0]*=f;
      if(!isnan(res.voi[ri].cl2[0])) res.voi[ri].cl2[0]*=f;
      if(!isnan(res.voi[ri].sd[0])) res.voi[ri].sd[0]*=f;
      if(per_dl) {
        res.voi[ri].parameter[0]*=100.;
        if(!isnan(res.voi[ri].cl1[0])) res.voi[ri].cl1[0]*=100.;
        if(!isnan(res.voi[ri].cl2[0])) res.voi[ri].cl2[0]*=100.;
        if(!isnan(res.voi[ri].sd[0])) res.voi[ri].sd[0]*=100.;
      }
      F=res.voi[ri].parameter[0]; // for printing
      /* pWater as such */
      pWater=res.voi[ri].parameter[1];
      /* Va */
      if(per_dl) {
        res.voi[ri].parameter[2]*=100.;
        if(!isnan(res.voi[ri].cl1[2])) res.voi[ri].cl1[2]*=100.;
        if(!isnan(res.voi[ri].cl2[2])) res.voi[ri].cl2[2]*=100.;
        if(!isnan(res.voi[ri].sd[2])) res.voi[ri].sd[2]*=100.;
      }
      Va=res.voi[ri].parameter[2];
      /* time delay as such */
      delayT=res.voi[ri].parameter[3];
      /* Print fit results */
      if(verbose>4)
        fprintf(stdout, "%-20s  %6.2f %6.3f %6.2f %6.2f %8.2e\n",
                ttac.voi[ri].name, F, pWater, Va, delayT, wss);
    }
 
    /* done with this region */
    if(ttac.voiNr>2 && verbose==1) {fprintf(stdout, "."); fflush(stdout);}
 
  } /* Next region */
  if(verbose>0) {fprintf(stdout, "\n"); fflush(stdout);}
 
  /*
   *  Print results on screen
   */
  if(verbose>0) {resPrint(&res); fprintf(stdout, "\n");}
 
 
  /*
   *  Save results
   */
  if(verbose>1) printf("saving results\n");
  if(resWrite(&res, resfile, verbose-3)!=0) {
    fprintf(stderr, "Error in writing '%s': %s\n", resfile, reserrmsg);
    dftEmpty(&btac); dftEmpty(&ttac); resEmpty(&res);
    return(11);
  }
  if(verbose>1) fprintf(stdout, "Model parameters written in %s\n", resfile);
 
 
  /*
   *  Saving and/or plotting of fitted TACs
   */
  if(svgfile[0] || fitfile[0]) {
 
    if(times_changed) dftSec2min(&ttac);
 
    /* Create a DFT containing fitted TACs */
    int ri, fi, ret;
    char tmp[64];
    DFT dft2;
    dftInit(&dft2); ret=dftdup(&ttac, &dft2);
    if(ret) {
      fprintf(stderr, "Error: cannot save fitted curves.\n");
      dftEmpty(&ttac); dftEmpty(&btac); resEmpty(&res);
      return(21);
    }
    for(ri=0; ri<ttac.voiNr; ri++) for(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, "Radiowater fit ");
      if(strlen(res.studynr)>0) strcat(tmp, res.studynr);
      ret=plot_fitrange_svg(&ttac, &dft2, tmp, 0.0, 1.02*dft2.x[fitframeNr-1],
                            0.0, nan(""), svgfile, verbose-8);
      if(ret) {
        fprintf(stderr, "Error (%d) in writing '%s'.\n", ret, svgfile);
        dftEmpty(&dft2); dftEmpty(&ttac); dftEmpty(&btac); resEmpty(&res);
        return(30+ret);
      }
      if(verbose>0) printf("Plots written in %s\n", svgfile);
    }
 
    /* Save fitted TACs */
    if(fitfile[0]) {
      if(verbose>1) printf("saving fitted curves\n");
      tpcProgramName(argv[0], 1, 0, tmp, 64);
      sprintf(dft2.comments, "# program := %s\n", tmp);    
      if(dftWrite(&dft2, fitfile)) {
        fprintf(stderr, "Error in writing '%s': %s\n", fitfile, dfterrmsg);
        dftEmpty(&dft2); dftEmpty(&ttac); dftEmpty(&btac); resEmpty(&res);
        return(22);
      }
      if(verbose>0) printf("Fitted TACs written in %s\n", fitfile);
    }
 
    dftEmpty(&dft2);
  }
 
  /*
   *  Save residual curves, if required
   */
  if(residfile[0]) {
    if(verbose>1) printf("saving weighted residuals\n");
    int ri, fi;
    char tmp[64];
    for(ri=0; ri<ttac.voiNr; ri++) for(fi=0; fi<fitframeNr; fi++)
      ttac.voi[ri].y[fi]=ttac.voi[ri].y3[fi];
    tpcProgramName(argv[0], 1, 0, tmp, 64);
    sprintf(ttac.comments, "# Residual curves from %s\n", tmp);
    ttac.frameNr=fitframeNr;
    if(dftWrite(&ttac, residfile)) {
      fprintf(stderr, "Error in writing '%s': %s\n", residfile, dfterrmsg);
      dftEmpty(&ttac); dftEmpty(&btac); resEmpty(&res);
      return(13);
    }
    if(verbose>0) fprintf(stdout, "Residual curves written in %s\n", residfile);
  }
 
 
  /*
   *  Calculate and save time delay corrected input curve, if required
   */
  if(inputfile[0]) {
    int ri, fi;
    double *darr, time_diff;
    char tmp[64];
    if(verbose>1) printf("saving time delay corrected input curve\n");
    /* Calculate the median of delay times */
    darr=(double*)calloc(res.voiNr, sizeof(double));
    for(ri=0; ri<res.voiNr; ri++) {
      darr[ri]=res.voi[ri].parameter[3];
      if(verbose>4) printf("darr[%d]=%g\n", ri, darr[ri]);
    }
    time_diff=dmedian(darr, res.voiNr);
    free(darr);
    if(verbose>1) printf("time_diff := %g\n", time_diff);
    /* Change the sample times */
    for(fi=0; fi<btac.frameNr; fi++) {
      btac.x[fi]+=time_diff; 
      btac.x1[fi]+=time_diff; btac.x2[fi]+=time_diff;
    }
    /* 'remove' extra columns */
    btac.voiNr=1;
    /* Replace old comments, because those may contain wrong units */
    dftSetComments(&btac);
    sprintf(tmp, "# time_delay := %g [s]\n", time_diff);
    strcat(btac.comments, tmp);
    /* Write the file */
    if(dftWrite(&btac, inputfile)) {
      fprintf(stderr, "Error in writing '%s': %s\n", inputfile, dfterrmsg);
      dftEmpty(&ttac); dftEmpty(&btac); resEmpty(&res);
      return(15);
    }
    if(verbose>0) fprintf(stdout, "Input TAC written in %s\n", inputfile);
  }
 
 
  /* Free memory */
  dftEmpty(&ttac); dftEmpty(&btac); resEmpty(&res);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
/*
 *  Function to be minimized
 */
double waterFunc(int parNr, double *p, void *fdata)
{
  int fi, ret;
  double wss, d, k1, pWater, k2, Va, delayT;
  double pa[MAX_PARAMETERS], penalty=1.0;
 
 
  /* Check parameters against the constraints */
  ret=modelCheckParameters(parNr, pmin, pmax, p, pa, &penalty);
  if(fdata) {}
 
  /* Get parameters */
  k1=pa[0]; //k1/=(1.0-pa[2]);
  pWater=pa[1]; Va=pa[2]; delayT=pa[3];
  if(pWater>0.0) k2=k1/pWater; else k2=0.; //k2-=lambda;
  //printf("f=%g pWater=%g Va=%g\n", 6000.*k1, pWater, 100.*Va); 
 
  /* Simulate the tissue curve (k1 and k2 per perfusable tissue) */
  ret=simC3vs(btac.x, ca, ca, btac.frameNr, k1, k2, 0., 0., 0., 0.,
              0.0, Va, 1.0,
              cs, NULL, NULL, NULL, NULL, NULL);
  if(ret) {
    printf("error %d in simulation\n", ret);
    return(nan(""));
  }
 
  /* Correct simulated tissue curve for delay */
  if(fabs(delayT)>1.0E-06) {
    double dx[btac.frameNr];
    for(fi=0; fi<btac.frameNr; fi++) dx[fi]=btac.x[fi]+delayT;
    ret=interpolate(dx, cs, btac.frameNr, btac.x,
                  cs2, NULL, NULL, btac.frameNr);
    if(ret) {
      printf("error %d in interpolation\n", ret);
      return(nan(""));
    }
    for(fi=0; fi<btac.frameNr; fi++) cs[fi]=cs2[fi];
  }
  /* Remove decay correction */
  for(fi=0; fi<btac.frameNr; fi++) cs[fi]*=simdc[fi];
  /* Interpolate simulated TAC to PET frames */
  if(ttac.timetype==DFT_TIME_STARTEND)
    ret=interpolate4pet(btac.x, cs, btac.frameNr, ttac.x1, ttac.x2,
                        csim, NULL, NULL, fitframeNr);
  else
    ret=interpolate(btac.x, cs, btac.frameNr, ttac.x,
                        csim, NULL, NULL, fitframeNr);
  if(ret) {
    printf("error %d in interpolation\n", ret);
    return(nan(""));
  }
  /* Correct simulated PET curve for physical decay */
  for(fi=0; fi<fitframeNr; fi++) csim[fi]*=petdc[fi];
 
  /* Calculate wss */
  for(fi=0, wss=0.; fi<fitframeNr; fi++) {
    /* residual */
    d=cpet[fi]-csim[fi];
    /* error */
    if(weight[fi]>0.0) wss+=weight[fi]*d*d;
    /* weighted residual */
    if(weight[fi]>0.0) wres[fi]=sqrt(weight[fi])*d; else wres[fi]=0.0;
  }
  wss_wo_penalty=wss;
  wss*=penalty;
 
  //printf("f=%g pWater=%g Va=%g\n", 6000.*k1, pWater, 100.*Va); 
 
  return(wss);
}
/*****************************************************************************/
 
/*****************************************************************************/