lhsol.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 "tpctacmod.h"
#include "tpclinopt.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Fitting of full or reduced compartmental model to plasma and tissue",
  "time-activity curves (PTAC and TTAC) to estimate the model parameters.",
  " ",
  "  ____        ____        ____        ____                             ",
  " | Cp |--K1->| C1 |--k3->| C2 |--k5->| C3 |  compartments in series (s)",
  " |____|<-k2--|____|<-k4--|____|<-k6--|____|                            ",
  " ",
  "              ____        ____  ",
  "  ____       |    |--k3->| C2 |            compartments in parallel (p)",
  " |    |--K1->|    |<-k4--|____| ",
  " | Cp |      | C1 |       ____  ",
  " |____|<-k2--|    |--k5->| C3 | ",
  "             |____|<-k6--|____| ",
  " ",
  "Compartmental models are transformed into general linear least squares",
  "functions (1, 2, 3, 4), which are solved using Lawson-Hanson linear",
  "least-squares algorithms (5). Note that rate constants and macroparameters",
  "are represented per volume (as measured by PET) including vascular volume.", 
  " ",
  "Usage: @P [options] PTAC TTAC fittime results",
  " ",
  "Options:",
  " -model=<k1 | k2 | k3 | k4 | k5s | k6s | k5p | k6p>",
  "     representing the following compartmental model settings:",
  "     k1 (for assuming k2=k3=k4=k5=k6=0)",
  "     k2 (for assuming k3=k4=k5=k6=0)",
  "     k3 (for assuming k4=k5=k6=0)",
  "     k4 (for assuming k5=k6=0); default",
  "     k5s (for assuming k6=0 and compartments in series)",
  "     k6s (compartments in series)",
  "     k5p (for assuming k6=0 and compartments in parallel)",
//  "     k6p (compartments in parallel)",
//  "     For model 'k6p' most model parameters cannot be solved.",
  " -Vp=<ignored|fitted>",
  "     Vascular volume is ignored (default) or fitted; note that PTAC is",
  "     assumed to represent vascular blood curve.",
  " -w1 | -wf",
  "     Sample weights are set to 1 (-w1) or to frame lengths (-wf);",
  "     by default weights in TTAC file are used, if available.",
  " -mid",
  "     Use frame mid times even when start and end times are available.",
  " -svg=<Filename>",
  "     Fitted and measured TACs are plotted in specified SVG file.",
  " -fit=<Filename>",
  "     Fitted regional TTACs are written in specified file.",
  " -lp=<Filename>",
  "     Parameters of linear model are saved in specified file.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "References:",
  "1. Blomqvist G. On the construction of functional maps in positron emission",
  "   tomography. J Cereb Blood Flow Metab 1984;4:629-632.",
  "2. Gjedde A, Wong DF. Modeling neuroreceptor binding of radioligands",
  "   in vivo. In: Quantitative imaging: neuroreceptors, neurotransmitters,",
  "   and enzymes. (Eds. Frost JJ, Wagner HM Jr). Raven Press, 1990, 51-79.",
  "3. Oikonen V. Multilinear solution for 4-compartment model:",
  "   I. Tissue compartments in series.",
  "   https://www.turkupetcentre.net/reports/tpcmod0023.pdf",
  "4. Oikonen V. Multilinear solution for 4-compartment model:",
  "   II. Two parallel tissue compartments.",
  "   https://www.turkupetcentre.net/reports/tpcmod0024.pdf",
  "5. Lawson CL & Hanson RJ. Solving least squares problems.",
  "   Prentice-Hall, 1974.",
  " ",
  "See also: fitk4, fitk5, patlak, logan, imglhdv, fitdelay, taccbv",
  " ",
  "Keywords: TAC, modelling, compartmental model, LLSQ",
  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;
/*****************************************************************************/
 
/*****************************************************************************/
#define MAX_LLSQ_N 7
enum {MODEL_UNKNOWN, MODEL_K1, MODEL_K2, MODEL_K3, MODEL_K4,
      MODEL_K5S, MODEL_K5P, MODEL_K6S, MODEL_K6P};
enum {VB_UNKNOWN, VB_IGNORED, VB_FITTED};
enum {METHOD_UNKNOWN, METHOD_NNLS, METHOD_BVLS};
static char *model_str[] = {
  "unknown",
  "K1", "K1-k2", "K1-k3", "K1-k4", "K1-k5", "K1-k5 parallel",
  "K1-k6", "K1-k6 parallel", 0};
static char *vb_model_str[] = {"unknown", "ignored", "fitted", 0};
static char *method_str[] = {"unknown", "NNLS", "BVLS", 0};
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int     ai, help=0, version=0, verbose=1;
  char    ptacfile[FILENAME_MAX], ttacfile[FILENAME_MAX], resfile[FILENAME_MAX],
          fitfile[FILENAME_MAX], svgfile[FILENAME_MAX], lpfile[FILENAME_MAX];
  int     weights=0; // 0=default, 1=no weighting, 2=frequency
  int     mid=0; // 0=default, 1=use frame times even if start and end times exist.
  double  fitdur=nan("");
  int     method=METHOD_NNLS;
  int     model=MODEL_K4; //MODEL_UNKNOWN;
  int     vb_model=VB_IGNORED;
  int     ret;
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  ptacfile[0]=ttacfile[0]=resfile[0]=fitfile[0]=svgfile[0]=lpfile[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(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, "LP=", 3)==0) {
      strlcpy(lpfile, cptr+3, FILENAME_MAX); if(strlen(lpfile)>0) continue;
    } else if(strcasecmp(cptr, "W1")==0) {
      weights=1; continue;
    } else if(strcasecmp(cptr, "WF")==0) {
      weights=2; continue;
    } else if(strcasecmp(cptr, "MID")==0) {
      mid=1; continue;
    } else if(strcasecmp(cptr, "NNLS")==0) {
      method=METHOD_NNLS; continue;
    } else if(strcasecmp(cptr, "BVLS")==0) {
      method=METHOD_BVLS; continue;
    } else if(strncasecmp(cptr, "VB=", 3)==0 || 
              strncasecmp(cptr, "VP=", 3)==0 ||
              strncasecmp(cptr, "VA=", 3)==0) 
    {
      cptr+=3;
      if(strncasecmp(cptr, "FITTED", 1)==0) {vb_model=VB_FITTED; continue;}
      if(strncasecmp(cptr, "IGNORED", 1)==0) {vb_model=VB_IGNORED; continue;}
    } else if(strncasecmp(cptr, "MODEL=", 6)==0) {
      cptr+=6;
      if(strcasecmp(cptr, "K1")==0) {model=MODEL_K1; continue;}
      if(strcasecmp(cptr, "K2")==0) {model=MODEL_K2; continue;}
      if(strcasecmp(cptr, "K3")==0) {model=MODEL_K3; continue;}
      if(strcasecmp(cptr, "K4")==0) {model=MODEL_K4; continue;}
      if(strcasecmp(cptr, "K5S")==0) {model=MODEL_K5S; continue;}
      if(strcasecmp(cptr, "K5P")==0) {model=MODEL_K5P; continue;}
      if(strcasecmp(cptr, "K6S")==0) {model=MODEL_K6S; continue;}
      if(strcasecmp(cptr, "K6P")==0) {model=MODEL_K6P; continue;}
      fprintf(stderr, "Error: invalid model '%s'.\n", cptr);
      return(1);
    }
    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(ptacfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(ttacfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {
    if(atofCheck(argv[ai], &fitdur)) {
      fprintf(stderr, "Error: invalid fit time '%s'.\n", argv[ai]);
      return(1);
    }
    if(fitdur<=0.0) fitdur=1.0E+99;
    ai++;
  }
  if(ai<argc) strlcpy(resfile, 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(!resfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("ptacfile := %s\n", ptacfile);
    printf("ttacfile := %s\n", ttacfile);
    printf("resfile := %s\n", resfile);
    if(fitfile[0]) printf("fitfile := %s\n", fitfile);
    if(svgfile[0]) printf("svgfile := %s\n", svgfile);
    if(lpfile[0]) printf("lpfile := %s\n", lpfile);
    printf("model := %s\n", model_str[model]);
    printf("vb_model := %s\n", vb_model_str[vb_model]);
    printf("method := %s\n", method_str[method]);
    printf("required_fittime := %g min\n", fitdur);
    printf("weights := %d\n", weights);
    if(mid!=0) printf("mid := %d\n", mid);
  }
 
  /*
   *  Set model-dependent parameters
   */
  int llsq_n=0;
  switch(model) {
    case MODEL_K1:  llsq_n=1; break;
    case MODEL_K2:  llsq_n=2; break;
    case MODEL_K3:  llsq_n=3; break;
    case MODEL_K4:  llsq_n=4; break;
    case MODEL_K5S: llsq_n=5; break;
    case MODEL_K6S: llsq_n=6; break;
    case MODEL_K5P: llsq_n=5; break;
    case MODEL_K6P: llsq_n=6; break;
    default: exit(1);
  }
  if(vb_model==VB_FITTED) llsq_n++;
  if(verbose>2) printf("llsq_n := %d\n", llsq_n);
 
 
  /*
   *  Read tissue and input data
   */
  if(verbose>1) printf("reading tissue and input data\n");
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  TAC ptac, ttac0; tacInit(&ptac); tacInit(&ttac0);
  int fitSampleNr;
  ret=tacReadModelingData(ttacfile, ptacfile, NULL, NULL, &fitdur, 0,
                          &fitSampleNr, &ttac0, &ptac, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac0); tacFree(&ptac); return(2);
  }
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(ttac0.format));
    printf("tacNr := %d\n", ttac0.tacNr);
    printf("ttac.sampleNr := %d\n", ttac0.sampleNr);
    printf("ptac.sampleNr := %d\n", ptac.sampleNr);
    printf("fitSampleNr := %d\n", fitSampleNr);
    printf("xunit := %s\n", unitName(ttac0.tunit));
    printf("yunit := %s\n", unitName(ttac0.cunit));
    printf("fitdur := %g s\n", fitdur);
  }
  if(mid>0) ttac0.isframe=ptac.isframe=0;
  if(fitSampleNr<llsq_n || ptac.sampleNr<llsq_n) {
    fprintf(stderr, "Error: too few samples in specified fit duration.\n");
    tacFree(&ttac0); tacFree(&ptac); return(2);
  }
//int origSampleNr=ttac0.sampleNr;
  ttac0.sampleNr=fitSampleNr;
 
  /* Add data weights, if requested */
  if(weights==1) {
    ttac0.weighting=WEIGHTING_OFF; 
    for(int i=0; i<ttac0.sampleNr; i++) ttac0.w[i]=1.0;
  } else if(weights==2) {
    if(tacWByFreq(&ttac0, ISOTOPE_UNKNOWN, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&ttac0); tacFree(&ptac); return(2);
    }
  } else if(!tacIsWeighted(&ttac0)) {
    if(verbose>0) fprintf(stderr, "Warning: data is not weighted.\n");
  }
 
 
  /*
   *  Interpolate and integrate PTAC to TTAC times
   */
  if(verbose>1) printf("integrating PTAC\n");
  TAC input0; tacInit(&input0); // interpolated PTAC
  TAC input1; tacInit(&input1); // 1st integral 
  TAC input2; tacInit(&input2); // 2nd integral
  TAC input3; tacInit(&input3); // 3rd integral
  ret=tacInterpolate(&ptac, &ttac0, &input0, &input1, &input2, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac0); tacFree(&ptac); 
    tacFree(&input0); tacFree(&input1); tacFree(&input2);
    return(3);
  }
  ret=tacDuplicate(&input2, &input3);
  if(ret==TPCERROR_OK)
    ret=liIntegrate(input2.x, input2.c[0].y, input2.sampleNr, input3.c[0].y, 3, 0);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: cannot make 3rd integral of PTAC.\n");
    tacFree(&ttac0);  tacFree(&ptac); 
    tacFree(&input0); tacFree(&input1); tacFree(&input2); tacFree(&input3);
    return(3);
  }
  /* Original PTAC should not be needed any more */
  tacFree(&ptac);
 
 
  /* Integrate TTAC */
  if(verbose>1) printf("integrating TTAC\n");
  TAC ttac1; tacInit(&ttac1); // 1st integral 
  TAC ttac2; tacInit(&ttac2); // 2nd integral
  TAC ttac3; tacInit(&ttac3); // 3rd integral
  ret=tacInterpolate(&ttac0, &ttac0, NULL, &ttac1, &ttac2, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac0);  tacFree(&ttac1);  tacFree(&ttac2);
    tacFree(&input0); tacFree(&input1); tacFree(&input2); tacFree(&input3);
  }
  ret=tacDuplicate(&ttac2, &ttac3);
  if(ret==TPCERROR_OK) {
    for(int i=0; i<ttac2.tacNr; i++) {
      ret=liIntegrate(ttac2.x, ttac2.c[i].y, ttac2.sampleNr, ttac3.c[i].y, 3, 0);
      if(ret!=TPCERROR_OK) break;
    }
  }
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: cannot make 3rd integral of TTAC.\n");
    tacFree(&ttac0);  tacFree(&ttac1);  tacFree(&ttac2);  tacFree(&ttac3);
    tacFree(&input0); tacFree(&input1); tacFree(&input2); tacFree(&input3);
    return(3);
  }
 
 
  /*
   *  Prepare the room for LLSQ parameters
   */
  if(verbose>1) printf("initializing LLSQ parameter data\n");
  PAR lp; parInit(&lp);
  ret=parAllocateWithTAC(&lp, &ttac0, MAX_LLSQ_N, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac0);  tacFree(&ttac1);  tacFree(&ttac2);  tacFree(&ttac3);
    tacFree(&input0); tacFree(&input1); tacFree(&input2); tacFree(&input3);
    return(4);
  }
  /* Copy titles & file names */
  {
    int i;
    char buf[256];
    time_t t=time(NULL);
    /* set program name */
    tpcProgramName(argv[0], 1, 1, buf, 256);
    iftPut(&lp.h, "program", buf, 0, NULL);
    /* set file names */
    iftPut(&lp.h, "plasmafile", ptacfile, 0, NULL);
    iftPut(&lp.h, "datafile", ttacfile, 0, NULL);
    /* Fit method */
    iftPut(&lp.h, "fitmethod", method_str[method], 0, NULL);
    /* Model */
    iftPut(&lp.h, "model", model_str[model], 0, NULL);
    /* Set current time to results */
    iftPut(&lp.h, "analysis_time", ctime_r_int(&t, buf), 0, NULL);
    /* Set fit times for each TAC */
    for(i=0; i<lp.tacNr; i++) {
      lp.r[i].dataNr=fitSampleNr;
      lp.r[i].start=0.0;
      lp.r[i].end=fitdur;
      /* and nr of fitted parameters */
      lp.r[i].fitNr=llsq_n;
    }
    /* Set the parameter names and units */
    for(i=0; i<MAX_LLSQ_N; i++) {
      sprintf(lp.n[i].name, "P%d", 1+i);
    }
    lp.parNr=llsq_n;
  }
 
 
  if(method==METHOD_NNLS) {
 
    /*
     *  Allocate memory required by NNLS
     */
    if(verbose>1) printf("allocating memory for NNLS\n");
    int llsq_m=fitSampleNr;
    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(&ttac0);  tacFree(&ttac1);  tacFree(&ttac2);  tacFree(&ttac3);
      tacFree(&input0); tacFree(&input1); tacFree(&input2); tacFree(&input3);
      parFree(&lp);
      return(5);
    }
    double **llsq_a=(double**)malloc(llsq_n*sizeof(double*));
    if(llsq_a==NULL) {
      fprintf(stderr, "Error: cannot allocate memory for NNLS.\n");
      tacFree(&ttac0);  tacFree(&ttac1);  tacFree(&ttac2);  tacFree(&ttac3);
      tacFree(&input0); tacFree(&input1); tacFree(&input2); tacFree(&input3);
      parFree(&lp); free(llsq_mat);
      return(5);
    }
    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;
 
    /*
     *  Fit each regional TTAC
     */
    for(int ti=0; ti<ttac0.tacNr; ti++) {
 
      if(verbose>1 && ttac0.tacNr>1) {
        printf("Region %d %s\n", 1+ti, ttac0.c[ti].name); fflush(stdout);}
 
      /* Setup data matrix A and vector B */
      for(int mi=0; mi<llsq_m; mi++)
        llsq_b[mi]=ttac0.c[ti].y[mi];
      int n=llsq_n; if(vb_model==VB_FITTED) n--;
      for(int mi=0; mi<llsq_m; mi++) {
        if(n>0) llsq_mat[mi]=input1.c[0].y[mi];           // PTAC 1st integral
        if(n>1) llsq_mat[mi+llsq_m]=-ttac1.c[ti].y[mi];   // TTAC 1st integral
        if(n>2) llsq_mat[mi+2*llsq_m]=input2.c[0].y[mi];  // PTAC 2nd integral
        if(n>3) llsq_mat[mi+3*llsq_m]=-ttac2.c[ti].y[mi]; // TTAC 2nd integral
        if(n>4) llsq_mat[mi+4*llsq_m]=input3.c[0].y[mi];  // PTAC 3rd integral
        if(n>5) llsq_mat[mi+5*llsq_m]=-ttac3.c[ti].y[mi]; // TTAC 3rd integral
        if(vb_model==VB_FITTED) 
          llsq_mat[mi+(llsq_n-1)*llsq_m]=input0.c[0].y[mi];  // PTAC 
      }
      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(&ttac0)) nnlsWght(llsq_n, llsq_m, llsq_a, llsq_b, ttac0.w);
      /* Compute NNLS */
      if(verbose>3) printf("starting NNLS...\n");
      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", ttac0.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", ttac0.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");
      }
      for(int ni=0; ni<llsq_n; ni++) lp.r[ti].p[ni]=llsq_x[ni];
      lp.r[ti].wss=r2;
      lp.r[ti].dataNr=tacWSampleNr(&ttac0);
      lp.r[ti].fitNr=llsq_n;
 
      /* Compute fitted TAC (into ttac1 since it is otherwise not needed) */
      for(int mi=0; mi<llsq_m; mi++) {
        ttac1.c[ti].y[mi]=0.0;
        for(int ni=0; ni<llsq_n; ni++) ttac1.c[ti].y[mi]+=llsq_x[ni]*matbackup[mi+ni*llsq_m];
      }
 
    } // next TAC
 
    /* Free allocated memory */
    free(llsq_a); free(llsq_mat);
 
  } else if(method==METHOD_BVLS) {
 
    /*
     *  Allocate memory required by BVLS
     */
    if(verbose>1) printf("allocating memory for BVLS\n");
    int llsq_m=fitSampleNr;
    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(&ttac0);  tacFree(&ttac1);  tacFree(&ttac2);  tacFree(&ttac3);
      tacFree(&input0); tacFree(&input1); tacFree(&input2); tacFree(&input3);
      parFree(&lp);
      return(5);
    }
    double b[llsq_m], x[MAX_LLSQ_N], bl[MAX_LLSQ_N], bu[MAX_LLSQ_N], w[llsq_n], zz[llsq_m];
    double act[llsq_m*(llsq_n+2)], r2;
    int istate[llsq_n+1], iterNr;
    double *matbackup=llsq_mat+llsq_n*llsq_m;
 
    /*
     *  Fit each regional TTAC
     */
    for(int ti=0; ti<ttac0.tacNr; ti++) {
 
      if(verbose>1 && ttac0.tacNr>1) {
        printf("Region %d %s\n", 1+ti, ttac0.c[ti].name); fflush(stdout);}
 
      /* Setup data matrix A and vector B */
      for(int mi=0; mi<llsq_m; mi++)
        b[mi]=ttac0.c[ti].y[mi];
      int n=llsq_n; if(vb_model==VB_FITTED) n--;
      for(int mi=0; mi<llsq_m; mi++) {
        if(n>0) llsq_mat[mi]=input1.c[0].y[mi];           // PTAC 1st integral
        if(n>1) llsq_mat[mi+llsq_m]=-ttac1.c[ti].y[mi];   // TTAC 1st integral
        if(n>2) llsq_mat[mi+2*llsq_m]=input2.c[0].y[mi];  // PTAC 2nd integral
        if(n>3) llsq_mat[mi+3*llsq_m]=-ttac2.c[ti].y[mi]; // TTAC 2nd integral
        if(n>4) llsq_mat[mi+4*llsq_m]=input3.c[0].y[mi];  // PTAC 3rd integral
        if(n>5) llsq_mat[mi+5*llsq_m]=-ttac3.c[ti].y[mi]; // TTAC 3rd integral
        if(vb_model==VB_FITTED) 
          llsq_mat[mi+(llsq_n-1)*llsq_m]=input0.c[0].y[mi];  // PTAC 
      }
      if(verbose>5) {
        printf("Matrix A and vector B:\n");
        for(int mi=0; mi<llsq_m; mi++) {
          printf("%.2e", llsq_mat[mi]);
          for(int ni=1; ni<llsq_n; ni++) printf(", %.2e", llsq_mat[mi+ni*llsq_m]);
          printf("; %.3e\n", 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(&ttac0)) llsqWght(llsq_n, llsq_m, NULL, llsq_mat, b, ttac0.w);
      /* Set istate vector to indicate that all parameters are non-bound */
      istate[llsq_n]=0; for(int ni=0; ni<llsq_n; ni++) istate[ni]=1+ni;
      /* Set parameter limits */
      if(vb_model==VB_FITTED) {bl[llsq_n-1]=0.0; bu[llsq_n-1]=1.0;}
      switch(model) {
        case MODEL_K1:
          bl[0]=0.0; bu[0]=10.0;
        break;
        case MODEL_K2:
          bl[0]=0.0; bu[0]=10.0;
          bl[1]=0.0; bu[1]=10.0;
          break;
        case MODEL_K3:
          bl[0]=0.0; bu[0]=10.0;
          bl[1]=0.0; bu[1]=10.0;
          bl[2]=0.0; bu[2]=10.0;
          break;
        case MODEL_K4:
          bl[0]=0.0; bu[0]=10.0;
          bl[1]=0.0; bu[1]=10.0;
          bl[2]=0.0; bu[2]=2.0;
          bl[3]=0.0; bu[3]=2.0;
          break;
        case MODEL_K5S:
          bl[0]=0.0; bu[0]=10.0;
          bl[1]=0.0; bu[1]=10.0;
          bl[2]=0.0; bu[2]=2.0;
          bl[3]=0.0; bu[3]=2.0;
          bl[4]=0.0; bu[4]=1.0;
          break;
        case MODEL_K5P:
          bl[0]=0.0; bu[0]=10.0;
          bl[1]=0.0; bu[1]=10.0;
          bl[2]=0.0; bu[2]=2.0;
          bl[3]=0.0; bu[3]=2.0;
          bl[4]=0.0; bu[4]=1.0;
          break;
        case MODEL_K6S:
          bl[0]=0.0; bu[0]=10.0;
          bl[1]=0.0; bu[1]=10.0;
          bl[2]=0.0; bu[2]=2.0;
          bl[3]=0.0; bu[3]=2.0;
          bl[4]=0.0; bu[4]=1.0;
          bl[5]=0.0; bu[5]=0.2;
          break;
        case MODEL_K6P:
          bl[0]=0.0; bu[0]=10.0;
          bl[1]=0.0; bu[1]=10.0;
          bl[2]=0.0; bu[2]=2.0;
          bl[3]=0.0; bu[3]=2.0;
          bl[4]=0.0; bu[4]=1.0;
          bl[5]=0.0; bu[5]=0.2;
          break;
        default: exit(1);
      }
      /* Set max iterations */
      iterNr=3*llsq_n;
      /* Compute BVLS */
      if(verbose>3) printf("starting BVLS...\n");
      ret=bvls(0, llsq_m, llsq_n, llsq_mat, b, bl, bu, x, w, act, zz, istate, &iterNr, verbose-3);
      if(verbose>3) printf("  ... done.\n");
      r2=w[0];
      if(ret!=0) {
        if(ret==-1) fprintf(stderr, "Warning: BVLS iteration max exceeded for %s\n", ttac0.c[ti].name);
        else fprintf(stderr, "Warning: no BVLS solution for %s\n", ttac0.c[ti].name);
        for(int ni=0; ni<llsq_n; ni++) x[ni]=0.0;
        r2=0.0;
      }
      if(verbose>4) {
        printf("solution_vector: %d", istate[0]);
        for(int ni=1; ni<llsq_n; ni++) printf(", %d", istate[ni]);
        printf("\n");
      }
      for(int ni=0; ni<llsq_n; ni++) lp.r[ti].p[ni]=x[ni];
      lp.r[ti].wss=r2;
      lp.r[ti].dataNr=tacWSampleNr(&ttac0);
      lp.r[ti].fitNr=llsq_n;
 
      /* Compute fitted TAC (into ttac1 since it is otherwise not needed) */
      for(int mi=0; mi<llsq_m; mi++) {
        ttac1.c[ti].y[mi]=0.0;
        for(int ni=0; ni<llsq_n; ni++) ttac1.c[ti].y[mi]+=x[ni]*matbackup[mi+ni*llsq_m];
      }
 
    } // next TAC
 
    /* Free allocated memory */
    free(llsq_mat);
 
  } else {
 
    fprintf(stderr, "Error: selected method not available.");
    tacFree(&ttac0);  tacFree(&ttac1);  tacFree(&ttac2);  tacFree(&ttac3);
    tacFree(&input0); tacFree(&input1); tacFree(&input2); tacFree(&input3);
    parFree(&lp);
    return(1);
 
  }
 
 
  /*
   *  Print original LLSQ parameters on screen
   */
  if(verbose>1 && lp.tacNr<50) 
    parWrite(&lp, stdout, PAR_FORMAT_TSV_UK /*PAR_FORMAT_RES*/, 0, &status);
 
 
  /*
   *  Save original LLSQ results, if requested
   */
  if(lpfile[0]) {
    if(verbose>1) printf("writing %s\n", lpfile);
    lp.format=parFormatFromExtension(lpfile);
    if(verbose>2) printf("result file format := %s\n", parFormattxt(lp.format));
    if(lp.format==PAR_FORMAT_UNKNOWN) lp.format=PAR_FORMAT_TSV_UK;
    FILE *fp; fp=fopen(lpfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing parameter file.\n");
      ret=TPCERROR_FAIL;
    } else {
      ret=parWrite(&lp, 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(&ttac0);  tacFree(&ttac1);  tacFree(&ttac2);  tacFree(&ttac3);
      tacFree(&input0); tacFree(&input1); tacFree(&input2); tacFree(&input3);
      parFree(&lp);
      return(11);
    }
    if(verbose>0) printf("Results saved in %s.\n", lpfile);
  }
 
 
 
  /*
   *  Prepare the room for CM parameters.
   *  Solve CM parameters from LLSQ parameters.
   */
  if(verbose>1) printf("initializing CM parameter data\n");
  PAR cmpar; parInit(&cmpar);
  ret=parAllocateWithTAC(&cmpar, &ttac0, MAX_LLSQ_N+4, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac0);  tacFree(&ttac1);  tacFree(&ttac2);  tacFree(&ttac3);
    tacFree(&input0); tacFree(&input1); tacFree(&input2); tacFree(&input3);
    parFree(&lp);
    return(21);
  }
  /* Copy titles & file names */
  {
    int i;
    char buf[256];
    time_t t=time(NULL);
    /* set program name */
    tpcProgramName(argv[0], 1, 1, buf, 256);
    iftPut(&cmpar.h, "program", buf, 0, NULL);
    /* set file names */
    iftPut(&cmpar.h, "plasmafile", ptacfile, 0, NULL);
    iftPut(&cmpar.h, "datafile", ttacfile, 0, NULL);
    /* Fit method */
    iftPut(&cmpar.h, "fitmethod", method_str[method], 0, NULL);
    /* Model */
    iftPut(&cmpar.h, "model", model_str[model], 0, NULL);
    /* Set current time to results */
    iftPut(&cmpar.h, "analysis_time", ctime_r_int(&t, buf), 0, NULL);
    /* Set fit times for each TAC */
    for(i=0; i<cmpar.tacNr; i++) {
      cmpar.r[i].dataNr=fitSampleNr;
      cmpar.r[i].start=0.0;
      cmpar.r[i].end=fitdur;
      /* and nr of fitted parameters */
      cmpar.r[i].fitNr=llsq_n;
    }
    /* Set the parameter names and units */
    cmpar.parNr=llsq_n;
    double k1, k2, k3, k4, k5, k6, k1k2, k3k4, k5k6, Ki, Vt, Vp;
    const double llimit=1.0E-06; // set to zero or NaN if lower than this
    switch(model) {
      case MODEL_K1:
        i=0; strcpy(cmpar.n[i].name,"K1"); //cmpar.n[i].unit=UNIT_ML_PER_ML_MIN;
        i=1; strcpy(cmpar.n[i].name,"Vp"); cmpar.n[i].unit=UNIT_PERCENTAGE;
        for(int ti=0; ti<cmpar.tacNr; ti++) {
          if(vb_model==VB_FITTED) Vp=lp.r[ti].p[llsq_n-1]; else Vp=0.0;
          cmpar.r[ti].p[0]=k1=lp.r[ti].p[0];
          cmpar.r[ti].p[1]=100.*Vp;
        }
        break;
      case MODEL_K2:
        i=0; strcpy(cmpar.n[i].name,"K1"); //cmpar.n[i].unit=UNIT_ML_PER_ML_MIN;
        i=1; strcpy(cmpar.n[i].name,"k2"); //cmpar.n[i].unit=UNIT_PER_MIN;
        i=2; strcpy(cmpar.n[i].name,"K1/k2"); //cmpar.n[i].unit=UNIT_ML_PER_ML;
        i=3; strcpy(cmpar.n[i].name,"Vp"); cmpar.n[i].unit=UNIT_PERCENTAGE;
        for(int ti=0; ti<cmpar.tacNr; ti++) {
          if(vb_model==VB_FITTED) Vp=lp.r[ti].p[llsq_n-1]; else Vp=0.0;
          k1=lp.r[ti].p[0]-Vp*lp.r[ti].p[1];
          k2=lp.r[ti].p[1];
          if(k1<llimit) k1=k2=0.0;
          if(k2<llimit) k2=0.0;
          k1k2=k1/k2;
          if(isfinite(k1)) cmpar.r[ti].p[0]=k1;
          if(isfinite(k2)) cmpar.r[ti].p[1]=k2;
          if(isfinite(k1k2)) cmpar.r[ti].p[2]=k1k2;
          cmpar.r[ti].p[3]=100.*Vp;
        }
        cmpar.parNr+=1;
        break;
      case MODEL_K3:
        i=0; strcpy(cmpar.n[i].name,"K1"); //cmpar.n[i].unit=UNIT_ML_PER_ML_MIN;
        i=1; strcpy(cmpar.n[i].name,"k2"); //cmpar.n[i].unit=UNIT_PER_MIN;
        i=2; strcpy(cmpar.n[i].name,"k3"); //cmpar.n[i].unit=UNIT_PER_MIN;
        i=3; strcpy(cmpar.n[i].name,"K1/k2"); //cmpar.n[i].unit=UNIT_ML_PER_ML;
        i=4; strcpy(cmpar.n[i].name,"Ki"); //cmpar.n[i].unit=UNIT_ML_PER_ML_MIN;
        i=5; strcpy(cmpar.n[i].name,"Vp"); cmpar.n[i].unit=UNIT_PERCENTAGE;
        for(int ti=0; ti<cmpar.tacNr; ti++) {
          if(vb_model==VB_FITTED) Vp=lp.r[ti].p[llsq_n-1]; else Vp=0.;
          k1=lp.r[ti].p[0]-Vp*lp.r[ti].p[1];
          k2=lp.r[ti].p[1]-lp.r[ti].p[2]/k1;
          k3=lp.r[ti].p[1]-k2;
          if(k1<llimit) k1=k2=k3=0.0;
          if(k2<llimit) k2=k3=0.0;
          if(k3<llimit) k3=0.0;
          k1k2=k1/k2;
          Ki=k1*k3/(k2+k3);
          if(isfinite(k1)) cmpar.r[ti].p[0]=k1;
          if(isfinite(k2)) cmpar.r[ti].p[1]=k2;
          if(isfinite(k3)) cmpar.r[ti].p[2]=k3;
          if(isfinite(k1k2)) cmpar.r[ti].p[3]=k1k2;
          if(isfinite(Ki)) cmpar.r[ti].p[4]=Ki;
          cmpar.r[ti].p[5]=100.*Vp;
        }
        cmpar.parNr+=2;
        break;
      case MODEL_K4:
        i=0; strcpy(cmpar.n[i].name,"K1"); //cmpar.n[i].unit=UNIT_ML_PER_ML_MIN;
        i=1; strcpy(cmpar.n[i].name,"k2"); //cmpar.n[i].unit=UNIT_PER_MIN;
        i=2; strcpy(cmpar.n[i].name,"k3"); //cmpar.n[i].unit=UNIT_PER_MIN;
        i=3; strcpy(cmpar.n[i].name,"k4"); //cmpar.n[i].unit=UNIT_PER_MIN;
        i=4; strcpy(cmpar.n[i].name,"K1/k2"); //cmpar.n[i].unit=UNIT_ML_PER_ML;
        i=5; strcpy(cmpar.n[i].name,"k3/k4"); //cmpar.n[i].unit=UNIT_UNITLESS;
        i=6; strcpy(cmpar.n[i].name,"Vt"); //cmpar.n[i].unit=UNIT_ML_PER_ML;
        i=7; strcpy(cmpar.n[i].name,"Vp"); cmpar.n[i].unit=UNIT_PERCENTAGE;
        for(int ti=0; ti<cmpar.tacNr; ti++) {
          if(vb_model==VB_FITTED) Vp=lp.r[ti].p[llsq_n-1]; else Vp=0.0;
          k1=lp.r[ti].p[0]-Vp*lp.r[ti].p[1];
          k2=lp.r[ti].p[1]-(lp.r[ti].p[2]-Vp*lp.r[ti].p[3])/k1;
          k3=lp.r[ti].p[1]-k2-lp.r[ti].p[3]/k2;
          k4=lp.r[ti].p[1]-k2-k3;
          if(k1<llimit) k1=k2=k3=k4=0.0;
          if(k2<llimit) k2=k3=k4=0.0;
          if(k3<llimit) k3=k4=0.0;
          if(k4<llimit) k4=0.0;
          k1k2=k1/k2;
          k3k4=k3/k4;
          Vt=k1k2*(1.0+k3k4);
          if(isfinite(k1)) cmpar.r[ti].p[0]=k1;
          if(isfinite(k2)) cmpar.r[ti].p[1]=k2;
          if(isfinite(k3)) cmpar.r[ti].p[2]=k3;
          if(isfinite(k4)) cmpar.r[ti].p[3]=k4;
          if(isfinite(k1k2)) cmpar.r[ti].p[4]=k1k2;
          if(isfinite(k3k4)) cmpar.r[ti].p[5]=k3k4;
          if(isfinite(Vt)) cmpar.r[ti].p[6]=Vt;
          cmpar.r[ti].p[7]=100.*Vp;
        }
        cmpar.parNr+=3;
        break;
      case MODEL_K5S:
        i=0; strcpy(cmpar.n[i].name,"K1"); //cmpar.n[i].unit=UNIT_ML_PER_ML_MIN;
        i=1; strcpy(cmpar.n[i].name,"k2"); //cmpar.n[i].unit=UNIT_PER_MIN;
        i=2; strcpy(cmpar.n[i].name,"k3"); //cmpar.n[i].unit=UNIT_PER_MIN;
        i=3; strcpy(cmpar.n[i].name,"k4"); //cmpar.n[i].unit=UNIT_PER_MIN;
        i=4; strcpy(cmpar.n[i].name,"k5"); //cmpar.n[i].unit=UNIT_PER_MIN;
        i=5; strcpy(cmpar.n[i].name,"K1/k2"); //cmpar.n[i].unit=UNIT_ML_PER_ML;
        i=6; strcpy(cmpar.n[i].name,"k3/k4"); //cmpar.n[i].unit=UNIT_UNITLESS;
        i=7; strcpy(cmpar.n[i].name,"Ki"); //cmpar.n[i].unit=UNIT_ML_PER_ML_MIN;
        i=8; strcpy(cmpar.n[i].name,"Vp"); cmpar.n[i].unit=UNIT_PERCENTAGE;
        for(int ti=0; ti<cmpar.tacNr; ti++) {
          if(vb_model==VB_FITTED) Vp=lp.r[ti].p[llsq_n-1]; else Vp=0.0;
          k1=lp.r[ti].p[0]-Vp*lp.r[ti].p[1];
          k2=lp.r[ti].p[1]-(lp.r[ti].p[2]-Vp*lp.r[ti].p[3])/k1;
          k3=lp.r[ti].p[1]-k2-(lp.r[ti].p[3]-lp.r[ti].p[4]/k1)/k2;
          k4=lp.r[ti].p[1]-k2-k3-(lp.r[ti].p[4]/k1)/k3;
          k5=lp.r[ti].p[1]-k2-k3-k4;
          if(k1<llimit) k1=k2=k3=k4=k5=0.0;
          if(k2<llimit) k2=k3=k4=k5=0.0;
          if(k3<llimit) k3=k4=k5=0.0;
          if(k4<llimit) k4=k5=0.0;
          if(k5<llimit) k5=0.0;
          k1k2=k1/k2;
          k3k4=k3/k4;
          Ki=k1*k3*k5/(k2*k4+k2*k5+k3*k5);
          if(isfinite(k1)) cmpar.r[ti].p[0]=k1;
          if(isfinite(k2)) cmpar.r[ti].p[1]=k2;
          if(isfinite(k3)) cmpar.r[ti].p[2]=k3;
          if(isfinite(k4)) cmpar.r[ti].p[3]=k4;
          if(isfinite(k5)) cmpar.r[ti].p[4]=k5;
          if(isfinite(k1k2)) cmpar.r[ti].p[5]=k1k2;
          if(isfinite(k3k4)) cmpar.r[ti].p[6]=k3k4;
          if(isfinite(Ki)) cmpar.r[ti].p[7]=Ki;
          cmpar.r[ti].p[8]=100.*Vp;
        }
        cmpar.parNr+=3;
        break;
      case MODEL_K5P:
        i=0; strcpy(cmpar.n[i].name,"K1"); //cmpar.n[i].unit=UNIT_ML_PER_ML_MIN;
        i=1; strcpy(cmpar.n[i].name,"k2"); //cmpar.n[i].unit=UNIT_PER_MIN;
        i=2; strcpy(cmpar.n[i].name,"k3"); //cmpar.n[i].unit=UNIT_PER_MIN;
        i=3; strcpy(cmpar.n[i].name,"k4"); //cmpar.n[i].unit=UNIT_PER_MIN;
        i=4; strcpy(cmpar.n[i].name,"k5"); //cmpar.n[i].unit=UNIT_PER_MIN;
        i=5; strcpy(cmpar.n[i].name,"K1/k2"); //cmpar.n[i].unit=UNIT_ML_PER_ML;
        i=6; strcpy(cmpar.n[i].name,"k3/k4"); //cmpar.n[i].unit=UNIT_UNITLESS;
        i=7; strcpy(cmpar.n[i].name,"Ki"); //cmpar.n[i].unit=UNIT_ML_PER_ML_MIN;
        i=8; strcpy(cmpar.n[i].name,"Vp"); cmpar.n[i].unit=UNIT_PERCENTAGE;
        for(int ti=0; ti<cmpar.tacNr; ti++) {
          if(vb_model==VB_FITTED) Vp=lp.r[ti].p[llsq_n-1]; else Vp=0.0;
          k1=lp.r[ti].p[0]-Vp*lp.r[ti].p[1];
          k2=lp.r[ti].p[1]-(lp.r[ti].p[2]-Vp*lp.r[ti].p[3])/k1;
          k4=(lp.r[ti].p[3]-lp.r[ti].p[4]/k1)/k2;
          k5=lp.r[ti].p[4]/(k1*k4);
          k3=lp.r[ti].p[1]-k2-k4-k5;
          if(k1<llimit) k1=k2=k3=k4=k5=0.0;
          if(k2<llimit) k2=k3=k4=k5=0.0;
          if(k3<llimit) k3=k4=0.0;
          if(k4<llimit) k4=0.0;
          if(k5<llimit) k5=0.0;
          k1k2=k1/k2;
          k3k4=k3/k4;
          Ki=k1*k5/(k2+k5);
          if(isfinite(k1)) cmpar.r[ti].p[0]=k1;
          if(isfinite(k2)) cmpar.r[ti].p[1]=k2;
          if(isfinite(k3)) cmpar.r[ti].p[2]=k3;
          if(isfinite(k4)) cmpar.r[ti].p[3]=k4;
          if(isfinite(k5)) cmpar.r[ti].p[4]=k5;
          if(isfinite(k1k2)) cmpar.r[ti].p[5]=k1k2;
          if(isfinite(k3k4)) cmpar.r[ti].p[6]=k3k4;
          if(isfinite(Ki)) cmpar.r[ti].p[7]=Ki;
          cmpar.r[ti].p[8]=100.*Vp;
        }
        cmpar.parNr+=3;
        break;
      case MODEL_K6S:
        i=0; strcpy(cmpar.n[i].name,"K1"); //cmpar.n[i].unit=UNIT_ML_PER_ML_MIN;
        i=1; strcpy(cmpar.n[i].name,"k2"); //cmpar.n[i].unit=UNIT_PER_MIN;
        i=2; strcpy(cmpar.n[i].name,"k3"); //cmpar.n[i].unit=UNIT_PER_MIN;
        i=3; strcpy(cmpar.n[i].name,"k4"); //cmpar.n[i].unit=UNIT_PER_MIN;
        i=4; strcpy(cmpar.n[i].name,"k5"); //cmpar.n[i].unit=UNIT_PER_MIN;
        i=5; strcpy(cmpar.n[i].name,"k6"); //cmpar.n[i].unit=UNIT_PER_MIN;
        i=6; strcpy(cmpar.n[i].name,"K1/k2"); //cmpar.n[i].unit=UNIT_ML_PER_ML;
        i=7; strcpy(cmpar.n[i].name,"k3/k4"); //cmpar.n[i].unit=UNIT_UNITLESS;
        i=8; strcpy(cmpar.n[i].name,"k5/k6"); //cmpar.n[i].unit=UNIT_UNITLESS;
        i=9; strcpy(cmpar.n[i].name,"Vt"); //cmpar.n[i].unit=UNIT_ML_PER_ML;
        i=10; strcpy(cmpar.n[i].name,"Vp"); cmpar.n[i].unit=UNIT_PERCENTAGE;
        for(int ti=0; ti<cmpar.tacNr; ti++) {
          if(vb_model==VB_FITTED) Vp=lp.r[ti].p[llsq_n-1]; else Vp=0.0;
          k1=lp.r[ti].p[0]-Vp*lp.r[ti].p[1];
          k2=lp.r[ti].p[1]-(lp.r[ti].p[2]-Vp*lp.r[ti].p[3])/k1;
          k3=lp.r[ti].p[1]-k2-(lp.r[ti].p[3]-(lp.r[ti].p[4]-Vp*lp.r[ti].p[5])/k1)/k2;
          k4=lp.r[ti].p[1]-k2-k3-((lp.r[ti].p[4]-Vp*lp.r[ti].p[5])/k1-lp.r[ti].p[5]/k2)/k3;
          k6=lp.r[ti].p[5]/(k2*k4);
          k5=lp.r[ti].p[1]-k2-k3-k4-k6;
          if(k1<llimit) k1=k2=k3=k4=k5=k6=0.0;
          if(k2<llimit) k2=k3=k4=k5=k6=0.0;
          if(k3<llimit) k3=k4=k5=k6=0.0;
          if(k4<llimit) k4=k5=k6=0.0;
          if(k5<llimit) k5=k6=0.0;
          if(k6<llimit) k6=0.0;
          k1k2=k1/k2;
          k3k4=k3/k4;
          k5k6=k5/k6;
          Vt=k1k2*(1.0+k3k4*(1.0+k5k6));
          if(isfinite(k1)) cmpar.r[ti].p[0]=k1;
          if(isfinite(k2)) cmpar.r[ti].p[1]=k2;
          if(isfinite(k3)) cmpar.r[ti].p[2]=k3;
          if(isfinite(k4)) cmpar.r[ti].p[3]=k4;
          if(isfinite(k5)) cmpar.r[ti].p[4]=k5;
          if(isfinite(k6)) cmpar.r[ti].p[5]=k6;
          if(isfinite(k1k2)) cmpar.r[ti].p[6]=k1k2;
          if(isfinite(k3k4)) cmpar.r[ti].p[7]=k3k4;
          if(isfinite(k5k6)) cmpar.r[ti].p[8]=k5k6;
          if(isfinite(Vt)) cmpar.r[ti].p[9]=Vt;
          cmpar.r[ti].p[10]=100.*Vp;
        }
        cmpar.parNr+=4;
        break;
      case MODEL_K6P:
        i=0; strcpy(cmpar.n[i].name,"K1"); //cmpar.n[i].unit=UNIT_ML_PER_ML_MIN;
        i=1; strcpy(cmpar.n[i].name,"k2"); //cmpar.n[i].unit=UNIT_PER_MIN;
        i=2; strcpy(cmpar.n[i].name,"k3"); //cmpar.n[i].unit=UNIT_PER_MIN;
        i=3; strcpy(cmpar.n[i].name,"k4"); //cmpar.n[i].unit=UNIT_PER_MIN;
        i=4; strcpy(cmpar.n[i].name,"k5"); //cmpar.n[i].unit=UNIT_PER_MIN;
        i=5; strcpy(cmpar.n[i].name,"k6"); //cmpar.n[i].unit=UNIT_PER_MIN;
        i=6; strcpy(cmpar.n[i].name,"K1/k2"); //cmpar.n[i].unit=UNIT_ML_PER_ML;
        i=7; strcpy(cmpar.n[i].name,"k3/k4"); //cmpar.n[i].unit=UNIT_UNITLESS;
        i=8; strcpy(cmpar.n[i].name,"k5/k6"); //cmpar.n[i].unit=UNIT_UNITLESS;
        i=9; strcpy(cmpar.n[i].name,"Vt"); //cmpar.n[i].unit=UNIT_ML_PER_ML;
        i=10; strcpy(cmpar.n[i].name,"Vp"); cmpar.n[i].unit=UNIT_PERCENTAGE;
        for(int ti=0; ti<cmpar.tacNr; ti++) {
          if(vb_model==VB_FITTED) Vp=lp.r[ti].p[llsq_n-1]; else Vp=0.0;
          k1=lp.r[ti].p[0]-Vp*lp.r[ti].p[1];
          k2=lp.r[ti].p[1]-(lp.r[ti].p[2]-Vp*lp.r[ti].p[3])/k1;
          if(k1<llimit) k1=k2=0.0;
          if(k2<llimit) k2=0.0;
          k3=k4=k5=k6=nan(""); // cannot be solved
          k1k2=k1/k2;
          k3k4=k5k6=nan(""); // cannot be solved
          Vt=cmpar.r[ti].p[4]/cmpar.r[ti].p[5]-Vp;
          if(isfinite(k1)) cmpar.r[ti].p[0]=k1;
          if(isfinite(k2)) cmpar.r[ti].p[1]=k2;
          if(isfinite(k1k2)) cmpar.r[ti].p[6]=k1k2;
          if(isfinite(Vt)) cmpar.r[ti].p[9]=Vt;
          cmpar.r[ti].p[10]=100.*Vp;
        }
        cmpar.parNr+=4;
        break;
      default: exit(1);
    }
    /* Copy R^2 etc */
    for(int ti=0; ti<cmpar.tacNr; ti++) {
      cmpar.r[ti].wss=lp.r[ti].wss;
      cmpar.r[ti].dataNr=lp.r[ti].dataNr;
      cmpar.r[ti].fitNr=lp.r[ti].fitNr;
    }
  }
 
  /*
   *  Print CM parameters on screen
   */
  if(verbose>0 && lp.tacNr<80) 
    parWrite(&cmpar, stdout, PAR_FORMAT_TSV_UK /*PAR_FORMAT_RES*/, 0, &status);
 
  /*
   *  Save CM results
   */
  {
    if(verbose>1) printf("writing %s\n", resfile);
    cmpar.format=parFormatFromExtension(resfile);
    if(verbose>2) printf("result file format := %s\n", parFormattxt(cmpar.format));
    if(cmpar.format==PAR_FORMAT_UNKNOWN) cmpar.format=PAR_FORMAT_TSV_UK;
    FILE *fp; fp=fopen(resfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing parameter file.\n");
      ret=TPCERROR_FAIL;
    } else {
      ret=parWrite(&cmpar, 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(&ttac0);  tacFree(&ttac1);  tacFree(&ttac2);  tacFree(&ttac3);
      tacFree(&input0); tacFree(&input1); tacFree(&input2); tacFree(&input3);
      parFree(&lp); parFree(&cmpar);
      return(22);
    }
    if(verbose>0) printf("Results saved in %s.\n", resfile);
  }
 
 
 
  /*
   *  SVG plot of fitted and original data
   */
  if(svgfile[0]) {
 
    if(verbose>1) printf("saving SVG plot\n");
    int i;
    char buf[128];
    sprintf(buf, "%s %s", method_str[method], model_str[model]);
    i=iftFindKey(&ttac0.h, "studynr", 0);
    if(i<0) i=iftFindKey(&ttac0.h, "study_number", 0);
    if(i>=0) {strcat(buf, ": "); strcat(buf, ttac0.h.item[i].value);}
    ret=tacPlotFitSVG(&ttac0, &ttac1, buf, 0.0, nan(""), 0.0, nan(""), svgfile, &status);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&ttac0);  tacFree(&ttac1);  tacFree(&ttac2);  tacFree(&ttac3);
      tacFree(&input0); tacFree(&input1); tacFree(&input2); tacFree(&input3);
      parFree(&lp); parFree(&cmpar);
      return(31);
    }
    if(verbose>0) printf("Plots written in %s.\n", svgfile);
  }
 
 
  /*
   *  Save fitted TTACs
   */
  if(fitfile[0]) {
    if(verbose>1) printf("writing %s\n", fitfile);
    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(&ttac1, fp, TAC_FORMAT_UNKNOWN, 1, &status);
      fclose(fp);
      if(ret!=TPCERROR_OK) fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    }
    if(ret!=TPCERROR_OK) {
      tacFree(&ttac0);  tacFree(&ttac1);  tacFree(&ttac2);  tacFree(&ttac3);
      tacFree(&input0); tacFree(&input1); tacFree(&input2); tacFree(&input3);
      parFree(&lp); parFree(&cmpar);
      return(32);
    }
    if(verbose>0) printf("fitted TACs saved in %s.\n", fitfile);
  }
 
 
  tacFree(&ttac0);  tacFree(&ttac1);  tacFree(&ttac2);  tacFree(&ttac3);
  tacFree(&input0); tacFree(&input1); tacFree(&input2); tacFree(&input3);
  parFree(&lp); parFree(&cmpar);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/