fitmtrap.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
/*****************************************************************************/
#include "tpcextensions.h"
#include "tpcift.h"
#include "tpctac.h"
#include "tpcpar.h"
#include "tpcli.h"
#include "tpcstatist.h"
#include "tpccm.h"
#include "tpctacmod.h"
#include "tpclinopt.h"
#include "tpcrand.h"
#include "tpcnlopt.h"
#include "tpcbfm.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Fitting of a compartmental model to regional myocardial PET data, where",
  "radiotracer is assumed to be trapped into cells at first pass,",
  "as chemical microspheres.",
  "The multilinear model is",
  "  Cpet(t) = Vb*Cb(t) + ((1-Vb)*K1+Vb*k2)*Integral[Cb(t)]",
  "            - k2*Integral[Cpet(t)]",
  " ",
  "BTAC from LV cavity is used as input function.",
  " ",
  "Usage: @P [Options] tacfile BTAC [parfile]",
  " ",
  "Options:",
  " -k2=<value> | -k2=<MTAC> | -k2=median",
  "     By default, full three-parameter model is fitted to all regional TACs",
  "     in tacfile are fitted independently. With these options, k2 can be",
  "     constrained to given value, or to value from a given large muscle",
  "     region, or to the median of k2 estimates of all regions.",
  " -w1 | -wf",
  "     By default, weights in data file are used, if available.",
  "     With these options all weights can be set to 1.0 (no weighting)",
  "     or based on frame durations.",
  " -svg=<Filename>",
  "     Fitted and measured TACs are plotted in specified SVG file.",
  " -fit=<Filename>",
  "     Fitted TACs are written in specified TAC file.",
  " -input=<Filename>",
  "     Input (1-HCT)*PTAC is written in specified file.",
  " -isvg=<Filename>",
  "     Input curves are plotted in specified SVG file.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: imgmtrap, b2ptrap, fitk2",
  " ",
  "Keywords: myocardium, perfusion",
  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     tacfile[FILENAME_MAX], parfile[FILENAME_MAX], svgfile[FILENAME_MAX];
  char     bname[MAX_TACNAME_LEN+1], mname[MAX_TACNAME_LEN+1];
  char     fitfile[FILENAME_MAX], icorfile[FILENAME_MAX], isvgfile[FILENAME_MAX];
  int      weights=0; // 0=default, 1=no weighting, 2=frequency
  double   fixedk2=nan("");
  /* Fix k2: 0=no; 1=fix to value given by user; 2=fix to median; 3=fix from large ROI */
  int      fixk2=0; 
 
 
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  tacfile[0]=parfile[0]=svgfile[0]=fitfile[0]=icorfile[0]=isvgfile[0]=(char)0;
  bname[0]=mname[0]=(char)0;
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    char *cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(!*cptr) continue;
    if(strcasecmp(cptr, "W1")==0) {
      weights=1; continue;
    } else if(strcasecmp(cptr, "WF")==0) {
      weights=2; continue;
    } else if(strncasecmp(cptr, "SVG=", 4)==0 && strlen(cptr)>4) {
      strlcpy(svgfile, cptr+4, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "ISVG=", 5)==0 && strlen(cptr)>4) {
      strlcpy(isvgfile, cptr+5, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "FIT=", 4)==0 && strlen(cptr)>4) {
      strlcpy(fitfile, cptr+4, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "INPUT=", 6)==0 && strlen(cptr)>4) {
      strlcpy(icorfile, cptr+6, FILENAME_MAX); continue;
    } else if(strcasecmp(cptr, "K2=MEDIAN")==0) {
      fixk2=2; fixedk2=nan(""); continue;
    } else if(strncasecmp(cptr, "K2=", 3)==0) {
      int ret=atofCheck(cptr+3, &fixedk2);
      if(ret==0) {fixk2=1;} 
      else {strlcpy(mname, cptr+3, MAX_TACNAME_LEN+1); fixk2=3;}
      continue;
    }
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    return(1);
  } else break;
 
  TPCSTATUS status; statusInit(&status);
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  status.verbose=verbose-3;
  
  /* 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(tacfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(bname, argv[ai++], MAX_TACNAME_LEN+1);
  if(ai<argc) strlcpy(parfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
  /* Did we get all the information that we need? */
  if(!tacfile[0] || !bname[0]) { // note that parameter file is optional
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("tacfile := %s\n", tacfile);
    printf("bname := %s\n", bname);
    printf("fixk2 := %d\n", fixk2);
    if(fixk2==1) printf("fixedk2 := %g\n", fixedk2);
    if(fixk2==3) printf("mname := %s\n", mname);
    if(parfile[0]) printf("parfile := %s\n", parfile);
    if(svgfile[0]) printf("svgfile := %s\n", svgfile);
    if(isvgfile[0]) printf("isvgfile := %s\n", isvgfile);
    if(fitfile[0]) printf("fitfile := %s\n", fitfile);
    if(icorfile[0]) printf("icorfile := %s\n", icorfile);
    printf("weights := %d\n", weights);
    fflush(stdout);
  }
 
 
  /*
   *  Read the data
   */
  if(verbose>1) printf("reading data\n");
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  TAC tac; tacInit(&tac);
  if(tacRead(&tac, tacfile, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); return(2);
  }
  if(verbose>1) {
    printf("fileformat := %s\n", tacFormattxt(tac.format));
    printf("tacNr := %d\n", tac.tacNr);
    printf("sampleNr := %d\n", tac.sampleNr);
    printf("xunit := %s\n", unitName(tac.tunit));
    printf("yunit := %s\n", unitName(tac.cunit));
    printf("isframe := %d\n", tac.isframe);
  }
  if(tac.tacNr<1) {
    fprintf(stderr, "Error: invalid file.\n");
    tacFree(&tac); return(2);
  } else if(tac.tacNr<2) {
    if(verbose>0) fprintf(stderr, "Error: file contains just one TAC.\n");
    tacFree(&tac); return(2);
  }
  if(tac.sampleNr<5) {
    fprintf(stderr, "Error: too few samples.\n");
    tacFree(&tac); return(2);
  }
  /* Check NaNs */
  if(tacNaNs(&tac)>0) {
    fprintf(stderr, "Error: data contains missing values.\n");
    tacFree(&tac); return(2);
  }
  /* Sort data by sample time */
  tacSortByTime(&tac, &status);
  /* Make sure that frame mid times are set, otherwise plotting may not work */
  tacSetX(&tac, NULL);
  /* Get x range */
  double xmin, xmax;
  if(tacXRange(&tac, &xmin, &xmax)!=0) {
    fprintf(stderr, "Error: invalid data sample times.\n");
    tacFree(&tac); return(2);
  }
  if(verbose>2) {
    printf("orig_xmin := %g\n", xmin);
    printf("orig_xmax := %g\n", xmax);
  }
  if(!(xmin>=0.0)) {
    fprintf(stderr, "Error: negative sample times.\n");
    tacFree(&tac); return(2);
  }
  /* Convert sample times into minutes */
  if(tac.tunit==UNIT_UNKNOWN) {
    if(xmax<10.0) tac.tunit=UNIT_MIN; else tac.tunit=UNIT_SEC;
  }
  if(tacXUnitConvert(&tac, UNIT_MIN, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); return(2);
  }
  if(tacXRange(&tac, &xmin, &xmax)!=0) {
    fprintf(stderr, "Error: invalid data sample times.\n");
    tacFree(&tac); return(2);
  }
  if(verbose>1) {
    printf("xmin := %g\n", xmin);
    printf("xmax := %g\n", xmax);
    printf("final_xunit := %s\n", unitName(tac.tunit));
  }
  /* Convert concentration units into kBq/mL, if possible */
  if(tac.cunit==UNIT_BQ_PER_ML) {
    tacYUnitConvert(&tac, UNIT_KBQ_PER_ML, NULL);
  }
 
 
 
  /* Select the LV BTAC */
  int ib=-1;
  {
    int n=tacSelectTACs(&tac, bname, 1, NULL);
    if(n<1) {
      fprintf(stderr, "Error: specified BTAC not found in %s.\n", tacfile);
      tacFree(&tac); return(3);
    } else if(n>1) {
      fprintf(stderr, "Error: %d TACs match '%s' in %s.\n", n, bname, tacfile);
      tacFree(&tac); return(3);
    }
    ib=tacFirstSelected(&tac);
    if(verbose>1) {
      printf("BTAC name := %s\n", tac.c[ib].name);
      printf("BTAC index := %d\n", ib);
    }
  }
 
  /* Select the Big muscle TTAC */
  int im=-1;
  if(fixk2==3) {
    int n=tacSelectTACs(&tac, mname, 1, NULL);
    if(n<1) {
      fprintf(stderr, "Error: specified MTAC not found in %s.\n", tacfile);
      tacFree(&tac); return(3);
    } else if(n>1) {
      fprintf(stderr, "Error: %d TACs match '%s' in %s.\n", n, mname, tacfile);
      tacFree(&tac); return(3);
    }
    im=tacFirstSelected(&tac);
    if(verbose>1) {
      printf("MTAC name := %s\n", tac.c[im].name);
      printf("MTAC index := %d\n", im);
    }
  }
 
  /* Add data weights, if requested */
  if(weights==1) {
    tac.weighting=WEIGHTING_OFF; 
    for(int i=0; i<tac.sampleNr; i++) tac.w[i]=1.0;
  } else if(weights==2) {
    if(tacWByFreq(&tac, ISOTOPE_UNKNOWN, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); return(3);
    }
  }
 
  /* Set place for TAC integrals / fitted TACs */
  TAC taci; tacInit(&taci);
 
 
  /*
   *  Prepare the room for results
   */
  if(verbose>1) printf("initializing parameter data\n");
  PAR par; parInit(&par);
  if(parAllocateWithTAC(&par, &tac, 3, &status) != TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); //tacFree(&taci);
    return(5);
  }
  /* Copy titles & file names */
  {
    char buf[256];
    time_t t=time(NULL);
    /* set program name */
    tpcProgramName(argv[0], 1, 1, buf, 256);
    iftPut(&par.h, "program", buf, 0, NULL);
    /* set file names */
    iftPut(&par.h, "datafile", tacfile, 0, NULL);
    /* Fit method */
    iftPut(&par.h, "fitmethod", "NNLS", 0, NULL);
    /* Model */
    iftPut(&par.h, "model", "trap", 0, NULL);
    /* Set current time to results */
    iftPut(&par.h, "analysis_time", ctime_r_int(&t, buf), 0, NULL);
    /* Fixed k2 */
    if(fixk2==1) iftPutDouble(&par.h, "k2", fixedk2, 0, NULL);
    else if(fixk2==2) iftPut(&par.h, "k2-constraint", "median", 0, NULL);
    else if(fixk2==3) iftPut(&par.h, "k2-constraint", tac.c[im].name, 0, NULL);
    /* Set fit times for each TAC */
    for(int i=0; i<par.tacNr; i++) {
      par.r[i].dataNr=tac.sampleNr;
      par.r[i].start=xmin;
      par.r[i].end=xmax;
      /* and nr of fitted parameters */
      par.r[i].fitNr=3; // May be set to 2 before saving if k2 is constrained
    }
    /* Set the parameter names and units */
    par.parNr=3; // May be set to 2 before saving if k2 is constrained
    strcpy(par.n[0].name, "Vb");
    strcpy(par.n[1].name, "K1");
    strcpy(par.n[2].name, "k2");
  }
  /* Delete BTAC from results before saving into file */
 
 
  /*
   *  Full 3-parameter model fit to all TACs if k2 is not fixed
   *  or to big muscle if that was asked for
   */
  if(fixk2!=1) {
    if(verbose>1) printf("full model fitting\n");
 
    /* TAC integrals */
    if(tacIntegrate(&tac, &taci, &status) != TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); return(5);
    }
 
    /* Allocate matrices for NNLS */
    int llsq_m=tac.sampleNr;
    int llsq_n=3;
    double *llsq_mat=(double*)malloc((2*llsq_n*llsq_m)*sizeof(double));
    if(llsq_mat==NULL) {
      fprintf(stderr, "Error: cannot allocate memory for NNLS.\n");
      tacFree(&tac);  tacFree(&taci); parFree(&par);
      return(6);
    }
    double **llsq_a=(double**)malloc(llsq_n*sizeof(double*));
    if(llsq_a==NULL) {
      tacFree(&tac);  tacFree(&taci); parFree(&par);
      return(6);
    }
    for(int ni=0; ni<llsq_n; ni++) llsq_a[ni]=llsq_mat+ni*llsq_m;
    double r2, llsq_b[llsq_m], llsq_x[llsq_n], llsq_wp[llsq_n], llsq_zz[llsq_m];
    int indexp[llsq_n];
    double *matbackup=llsq_mat+llsq_n*llsq_m;
 
    for(int ti=0; ti<tac.tacNr; ti++) {
 
      /* BTAC obviously not fitted */
      if(ti==ib) continue;
      /* If requested, none other than big muscle fitted */
      if(fixk2==3 && ti!=im) continue;
      if(verbose>1) {printf("Region %d %s\n", 1+ti, tac.c[ti].name); fflush(stdout);}
 
      /* Setup data matrix A and vector B */
      for(int mi=0; mi<llsq_m; mi++)
        llsq_b[mi]=tac.c[ti].y[mi];              // TTAC
      for(int mi=0; mi<llsq_m; mi++) {
        llsq_mat[mi]=tac.c[ib].y[mi];            // BTAC
        llsq_mat[mi+llsq_m]=taci.c[ib].y[mi];    // BTAC integral
        llsq_mat[mi+2*llsq_m]=-taci.c[ti].y[mi]; // TTAC integral
      }
      if(verbose>5) {
        printf("Matrix A and vector B:\n");
        for(int mi=0; mi<llsq_m; mi++) {
          printf("%.2e", llsq_a[0][mi]);
          for(int ni=1; ni<llsq_n; ni++) printf(", %.2e", llsq_a[ni][mi]);
          printf("; %.3e\n", llsq_b[mi]);
        }
      }
      /* Make a copy of A matrix for later use */
      for(int i=0; i<llsq_n*llsq_m; i++) matbackup[i]=llsq_mat[i];
 
      /* Apply data weights */
      if(tacIsWeighted(&tac)) nnlsWght(llsq_n, llsq_m, llsq_a, llsq_b, tac.w);
 
      /* Compute NNLS */
      if(verbose>3) printf("starting NNLS...\n");
      int ret=nnls(llsq_a, llsq_m, llsq_n, llsq_b, llsq_x, &r2, llsq_wp, llsq_zz, indexp);
      if(verbose>3) printf("  ... done.\n");
      if(ret>1) {
        fprintf(stderr, "Warning: no NNLS solution for %s\n", tac.c[ti].name);
        for(int ni=0; ni<llsq_n; ni++) llsq_x[ni]=0.0;
        r2=0.0;
      } else if(ret==1) {
        fprintf(stderr, "Warning: NNLS iteration max exceeded for %s\n", tac.c[ti].name);
      }
      if(verbose>4) {
        printf("solution_vector: %g", llsq_wp[0]);
        for(int ni=1; ni<llsq_n; ni++) printf(", %g", llsq_wp[ni]);
        printf("\n");
        printf("x coefficient_vector: %g", llsq_x[0]);
        for(int ni=1; ni<llsq_n; ni++) printf(", %g", llsq_x[ni]);
        printf("\n");
      }
      for(int ni=0; ni<llsq_n; ni++) par.r[ti].p[ni]=llsq_x[ni];
      par.r[ti].wss=r2;
      par.r[ti].dataNr=tacWSampleNr(&tac);
      par.r[ti].fitNr=llsq_n;
 
      /* Compute fitted TAC (into taci since it is otherwise not needed) */
      for(int mi=0; mi<llsq_m; mi++) {
        taci.c[ti].y[mi]=0.0;
        for(int ni=0; ni<llsq_n; ni++) taci.c[ti].y[mi]+=llsq_x[ni]*matbackup[mi+ni*llsq_m];
      }
    } // next TAC
 
    /* Free allocated memory */
    free(llsq_a); free(llsq_mat);
 
    /* Print original LLSQ parameters on screen */
    if(verbose>3 && par.tacNr<50) 
      parWrite(&par, stdout, PAR_FORMAT_TSV_UK, 0, &status);
 
    /* Convert LLSQ parameters to CM model parameters */
    /* Calculate median of k2 values */
    double k2list[tac.tacNr];
    int k2nr=0;
    for(int ti=0; ti<par.tacNr; ti++) {
      if(ti==ib) continue;
      if(fixk2==3 && ti!=im) continue;
      double p1=par.r[ti].p[0];
      double p2=par.r[ti].p[1];
      double p3=par.r[ti].p[2];
      p2-=p1*p3; if(p1<0.999) p2/=(1.0-p1);
      if(p2>=0.0) par.r[ti].p[1]=p2; else par.r[ti].p[1]=0.0;
      if(par.r[ti].p[0]<0.80 && par.r[ti].p[1]>0.0) k2list[k2nr++]=p3;
    }
    if(k2nr==1) {
      fixedk2=k2list[0];
    } else if(k2nr>1) {
      fixedk2=statMedian(k2list, k2nr);
    }
 
  }
 
  /*
   *  Simple 2-parameter model fit to all TACs with fixed k2, if requested
   */
  if(fixk2!=0) {
 
    if(verbose>1) printf("model fitting with k2 fixed to %g\n", fixedk2);
    if(!(fixedk2>=0.0)) {
      fprintf(stderr, "Error: invalid k2.\n");
      tacFree(&tac); tacFree(&taci); parFree(&par);
      return(7);
    }
 
    /* TAC integrals */
    if(tacIntegrate(&tac, &taci, &status) != TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); return(5);
    }
 
    /* Allocate matrices for NNLS */
    int llsq_m=tac.sampleNr;
    int llsq_n=2;
    double *llsq_mat=(double*)malloc((2*llsq_n*llsq_m)*sizeof(double));
    if(llsq_mat==NULL) {
      fprintf(stderr, "Error: cannot allocate memory for NNLS.\n");
      tacFree(&tac);  tacFree(&taci); parFree(&par);
      return(6);
    }
    double **llsq_a=(double**)malloc(llsq_n*sizeof(double*));
    if(llsq_a==NULL) {
      tacFree(&tac);  tacFree(&taci); parFree(&par);
      return(6);
    }
    for(int ni=0; ni<llsq_n; ni++) llsq_a[ni]=llsq_mat+ni*llsq_m;
    double r2, llsq_b[llsq_m], llsq_x[llsq_n], llsq_wp[llsq_n], llsq_zz[llsq_m];
    int indexp[llsq_n];
    double *matbackup=llsq_mat+llsq_n*llsq_m;
 
    for(int ti=0; ti<tac.tacNr; ti++) {
 
      /* BTAC obviously not fitted */
      if(ti==ib) continue;
      if(verbose>1) {printf("Region %d %s\n", 1+ti, tac.c[ti].name); fflush(stdout);}
 
      /* Setup data matrix A and vector B */
      for(int mi=0; mi<llsq_m; mi++)
        llsq_b[mi]=tac.c[ti].y[mi]+fixedk2*taci.c[ti].y[mi];   // TTAC + k2 * TTAC integral
      for(int mi=0; mi<llsq_m; mi++) {
        llsq_mat[mi]=tac.c[ib].y[mi]+fixedk2*taci.c[ib].y[mi]; // BTAC + k2 * BTAC integral
        llsq_mat[mi+llsq_m]=taci.c[ib].y[mi];                  // BTAC integral
      }
      if(verbose>5) {
        printf("Matrix A and vector B:\n");
        for(int mi=0; mi<llsq_m; mi++) {
          printf("%.2e", llsq_a[0][mi]);
          for(int ni=1; ni<llsq_n; ni++) printf(", %.2e", llsq_a[ni][mi]);
          printf("; %.3e\n", llsq_b[mi]);
        }
      }
      /* Make a copy of A matrix for later use */
      for(int i=0; i<llsq_n*llsq_m; i++) matbackup[i]=llsq_mat[i];
 
      /* Apply data weights */
      if(tacIsWeighted(&tac)) nnlsWght(llsq_n, llsq_m, llsq_a, llsq_b, tac.w);
 
      /* Compute NNLS */
      if(verbose>3) printf("starting NNLS...\n");
      int ret=nnls(llsq_a, llsq_m, llsq_n, llsq_b, llsq_x, &r2, llsq_wp, llsq_zz, indexp);
      if(verbose>3) printf("  ... done.\n");
      if(ret>1) {
        fprintf(stderr, "Warning: no NNLS solution for %s\n", tac.c[ti].name);
        for(int ni=0; ni<llsq_n; ni++) llsq_x[ni]=0.0;
        r2=0.0;
      } else if(ret==1) {
        fprintf(stderr, "Warning: NNLS iteration max exceeded for %s\n", tac.c[ti].name);
      }
      if(verbose>4) {
        printf("solution_vector: %g", llsq_wp[0]);
        for(int ni=1; ni<llsq_n; ni++) printf(", %g", llsq_wp[ni]);
        printf("\n");
        printf("x coefficient_vector: %g", llsq_x[0]);
        for(int ni=1; ni<llsq_n; ni++) printf(", %g", llsq_x[ni]);
        printf("\n");
      }
      for(int ni=0; ni<llsq_n; ni++) par.r[ti].p[ni]=llsq_x[ni];
      par.r[ti].wss=r2;
      par.r[ti].dataNr=tacWSampleNr(&tac);
      par.r[ti].fitNr=llsq_n;
 
      /* Compute fitted TAC (into taci since it is otherwise not needed) */
      /* Subtract k2 * TTAC integral so that fitted TACs are comparable to measured TACs */
      for(int mi=0; mi<llsq_m; mi++) {
        taci.c[ti].y[mi]*=-fixedk2;
        for(int ni=0; ni<llsq_n; ni++) taci.c[ti].y[mi]+=llsq_x[ni]*matbackup[mi+ni*llsq_m];
      }
 
    } // next TAC
 
    /* Free allocated memory */
    free(llsq_a); free(llsq_mat);
 
    /* Print original LLSQ parameters on screen */
    if(verbose>3 && par.tacNr<50) 
      parWrite(&par, stdout, PAR_FORMAT_TSV_UK, 0, &status);
 
    /* Convert LLSQ parameters to CM model parameters */
    for(int ti=0; ti<par.tacNr; ti++) {
      if(ti==ib) continue;
      double p1=par.r[ti].p[0];
      double p2=par.r[ti].p[1];
      if(p1<0.999) p2/=(1.0-p1);
      if(p2>=0.0) par.r[ti].p[1]=p2; else par.r[ti].p[1]=0.0;
      par.r[ti].p[2]=fixedk2;
    }
  }
 
 
  /*
   *  Calculate input curves for plotting/saving:
   *  Cb, (1-HCT)*Cp, and HCT*Crbc
   */
  if(icorfile[0] || isvgfile[0]) {
    TAC inp; tacInit(&inp);
    if(tacExtract(&tac, &inp, ib)!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot plot input curves.\n");
      tacFree(&tac);  tacFree(&taci); parFree(&par);
      return(101);
    }
    if(tacAllocateMore(&inp, 2)!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot allocate memory.\n");
      tacFree(&tac);  tacFree(&taci); parFree(&par);
      return(102);
    }
    inp.tacNr=3;
    strcpy(inp.c[0].name, "Blood");
    strcpy(inp.c[1].name, "Plasma");
    strcpy(inp.c[2].name, "RBC");
    /* Simulate (1-HCT)*Cp curve */
    double *cb=inp.c[0].y;
    double *cp=inp.c[1].y;
    double cpi=0.0;
    double t_last=0.0, cp_last=0.0, cpi_last=0.0;
    for(int i=0; i<inp.sampleNr; i++) {
      /* delta time / 2 */
      double dt2=0.5*(inp.x[i]-t_last);
      if(dt2>0.0) {
        /* plasma */
        cp[i] = (cb[i] - fixedk2*(cpi_last + dt2*cp_last)) / (1.0 + dt2*fixedk2);
        /* integral */
        cpi+=(cp[i]+cp_last)*dt2;
      } else {cp[i]=cp_last;}
      t_last=inp.x[i];
      cp_last=cp[i];
      cpi_last=cpi;
    }
    /* Calculate HCT*Crbc */
    for(int i=0; i<inp.sampleNr; i++) {
      inp.c[2].y[i]=inp.c[0].y[i] - inp.c[1].y[i];
    }
    if(isvgfile[0]) {
      if(verbose>1) printf("plotting input TAC\n");
      if(tacPlotLineSVG(&inp, "", 0.0, nan(""), nan(""), nan(""), isvgfile, NULL)!=TPCERROR_OK) {
        fprintf(stderr, "Error: cannot plot input curves.\n");
        tacFree(&tac);  tacFree(&taci); parFree(&par); tacFree(&inp);
        return(104);
      }
      if(verbose>0) printf("Input curves plotted in %s.\n", isvgfile);
    }
    if(icorfile[0]) {
      if(verbose>1) printf("saving (1-HCT)*Cp TAC\n");
      tacSwapTACCs(&inp, 0, 1);
      inp.tacNr=1;
      int ret=TPCERROR_OK;
      FILE *fp; fp=fopen(icorfile, "w");
      if(fp==NULL) {
        fprintf(stderr, "Error: cannot open file for writing '%s'.\n", icorfile);
        ret=TPCERROR_FAIL;
      } else {
        ret=tacWrite(&inp, fp, TAC_FORMAT_PMOD, 1, &status);
        fclose(fp);
        if(ret!=TPCERROR_OK) fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      }
      if(ret!=TPCERROR_OK) {
        tacFree(&tac);  tacFree(&taci); parFree(&par); tacFree(&inp);
        return(105);
      }
      if(verbose>0) printf("Corrected input saved in %s.\n", icorfile);
      inp.tacNr=3;
    }
    tacFree(&inp);
  }
 
 
  /*
   *  Save and/or plot fitted TACs, if requested
   */
  if(svgfile[0] || fitfile[0]) {
 
    /* Since fitted TACs are now in place of integral curves, but BTAC obviously was not fitted,
       copy the original BTAC in there, so that plots won't be messed up */
    for(int i=0; i<taci.sampleNr; i++) taci.c[ib].y[i]=tac.c[ib].y[i];
 
    /*
     *  SVG plot of fitted and original data
     */
    if(svgfile[0]) {
      if(verbose>1) printf("saving SVG plot\n");
      char buf[128];
      sprintf(buf, "NNLS trap");
      int i=iftFindKey(&tac.h, "studynr", 0);
      if(i<0) i=iftFindKey(&tac.h, "study_number", 0);
      if(i>=0) {strcat(buf, ": "); strcat(buf, tac.h.item[i].value);}
      if(tacPlotFitSVG(&tac, &taci, buf, 0.0, nan(""), 0.0, nan(""), svgfile, &status)!=TPCERROR_OK) {
        fprintf(stderr, "Error: %s\n", errorMsg(status.error));
        tacFree(&tac);  tacFree(&taci); parFree(&par);
        return(201);
      }
      if(verbose>0) printf("Plots written in %s.\n", svgfile);
    }
 
    /*
     *  Save fitted TTACs
     */
    if(fitfile[0]) {
      if(verbose>1) printf("writing %s\n", fitfile);
      int ret=TPCERROR_OK;
      FILE *fp; fp=fopen(fitfile, "w");
      if(fp==NULL) {
        fprintf(stderr, "Error: cannot open file for writing fitted TTACs.\n");
        ret=TPCERROR_FAIL;
      } else {
        ret=tacWrite(&taci, fp, TAC_FORMAT_PMOD, 1, &status);
        fclose(fp);
        if(ret!=TPCERROR_OK) fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      }
      if(ret!=TPCERROR_OK) {
        tacFree(&tac);  tacFree(&taci); parFree(&par);
        return(202);
      }
      if(verbose>0) printf("fitted TACs saved in %s.\n", fitfile);
    }
 
  }
 
 
  /* Remove BTAC from parameters */
  parDeleteTAC(&par, ib);
 
  /* If k2 was constrained to median or large muscle region, write it in header
     and remove it from parameters */
  if(fixk2==2 || fixk2==3) iftPutDouble(&par.h, "k2", fixedk2, 0, NULL);
  if(fixk2!=0) {
    par.parNr=2;
    for(int i=0; i<par.tacNr; i++) par.r[i].fitNr=2;
  }
 
  /* Show the parameters */
  if(verbose>0) parWrite(&par, stdout, PAR_FORMAT_TSV_UK, 1, NULL);
 
  /* 
   *  Save the parameters, if requested
   */
  if(parfile[0]) {
    par.format=parFormatFromExtension(parfile);
    if(verbose>2) printf("parameter file format := %s\n", parFormattxt(par.format));
    if(par.format==PAR_FORMAT_UNKNOWN) par.format=PAR_FORMAT_TSV_UK;
    /* Save file */
    if(verbose>1) printf("  saving %s\n", parfile);
    int ret=TPCERROR_OK;
    FILE *fp; fp=fopen(parfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing results.\n");
      ret=TPCERROR_FAIL;
    } else {
      ret=parWrite(&par, fp, PAR_FORMAT_UNKNOWN, 1, &status);
      fclose(fp);
      if(ret!=TPCERROR_OK) fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    }
    if(ret!=TPCERROR_OK) {
      tacFree(&tac);  tacFree(&taci); parFree(&par);
      return(21);
    }
    if(verbose>0) printf("parameters saved in %s\n", parfile);
  }
 
 
  parFree(&par);
  tacFree(&taci);
  tacFree(&tac);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/