llsqrk2.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[] = {
  "LLSQ fitting of reference tissue input model in case of one tissue",
  "compartment. Compartmental model is transformed into general linear least",
  "squares functions.",
  " ",
  "Radioactivity concentration in region(s)-of-interest is given in TTAC file;",
  "data must be corrected for physical decay.",
  "Radioactivity concentration in reference region (RTAC) is given in separate",
  "file, or as name or number of the reference region inside the TTAC file.",
  "Sample times must be in minutes in all data files, unless specified inside",
  "the files. TTAC file should include weights.",
  " ",
  "Usage: @P [options] TTAC RTAC results",
  " ",
  "Options:",
  " -end=<Fit end time (min)>",
  "     Use data from 0 to end time; by default, model is fitted to all frames.",
  " -w1 | -wf",
  "     Sample weights are set to 1 (-w1) or to frame lengths (-wf);",
  "     by default weights in TTAC file are used, if available.",
  " -model=<R1 | R2>",
  "     Use two-parameter model (R1, default), or three-parameter model (R2).",
/*
  " -k2r=<value>",
  "     k2r is fixed to given value, three-parameter model is set, and BVLS",
  "     is used instead of NNLS.",
  "     HIDDEN because does not yet work correctly but fixed R1*k2r instead.",
*/
  " -svg=<Filename>",
  "     Fitted and measured TACs are plotted in specified SVG file.",
  " -fit=<Filename>",
  "     Fitted regional TTACs are written 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. Lawson CL & Hanson RJ. Solving least squares problems.",
  "   Prentice-Hall, 1974.",
  " ",
  "See also: lhsol, bfmsrtm, 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 3
enum {METHOD_UNKNOWN, METHOD_NNLS, METHOD_BVLS};
static char *method_str[] = {"unknown", "NNLS", "BVLS", 0};
/*****************************************************************************/
enum {MODEL_UNKNOWN, MODEL_R1, MODEL_R1R2};
static char *model_str[] = {"unknown", "R1", "R1R2", 0};
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int     ai, help=0, version=0, verbose=1;
  char    ttacfile[FILENAME_MAX], rtacfile[FILENAME_MAX], resfile[FILENAME_MAX],
          fitfile[FILENAME_MAX], svgfile[FILENAME_MAX];
  int     weights=0; // 0=default, 1=no weighting, 2=frequency
  double  fitdur=nan("");
  double  fixed_k2r=nan("");
  int     method=METHOD_UNKNOWN;
  int     model=MODEL_UNKNOWN;
  int     ret;
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  ttacfile[0]=rtacfile[0]=resfile[0]=fitfile[0]=svgfile[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(strcasecmp(cptr, "W1")==0) {
      weights=1; continue;
    } else if(strcasecmp(cptr, "WF")==0) {
      weights=2; continue;
    } else if(strncasecmp(cptr, "END=", 4)==0) {
      if(atofCheck(cptr+4, &fitdur)==0) {
        if(fitdur<=0.0) fitdur=1.0E+99;
        continue;
      } else {
        fprintf(stderr, "Error: invalid fit time '%s'.\n", argv[ai]);
        return(1);
      }
    } else if(strncasecmp(cptr, "K2R=", 4)==0) {
      if(atofCheck(cptr+4, &fixed_k2r)==0 && fixed_k2r>0.0) continue;
      fprintf(stderr, "Error: invalid k2r '%s'.\n", argv[ai]); return(1);
    } else if(strcasecmp(cptr, "NNLS")==0) {
      method=METHOD_NNLS; continue;
    } else if(strcasecmp(cptr, "BVLS")==0) {
      method=METHOD_BVLS; continue;
    } else if(strncasecmp(cptr, "MODEL=", 6)==0) {
      cptr+=6;
      if(strcasecmp(cptr, "R1")  ==0) {model=MODEL_R1;   continue;}
      if(strcasecmp(cptr, "R2")  ==0) {model=MODEL_R1R2; continue;}
      if(strcasecmp(cptr, "R1R2")==0) {model=MODEL_R1R2; 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(ttacfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(rtacfile, argv[ai++], FILENAME_MAX);
  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);
  }
  /* Set default parameter values, if not set by user */
  if(isnan(fitdur)) fitdur=1.0E+99;
  if(method==METHOD_UNKNOWN) method=METHOD_NNLS;
  if(model==MODEL_UNKNOWN) model=MODEL_R1;
  /* If user fixed k2r, then model should include it */
  if(!isnan(fixed_k2r) && model==MODEL_R1) {
    model=MODEL_R1R2;
  }
  /* If user fixed k2r, then method should allow it */
  if(!isnan(fixed_k2r) && method==METHOD_NNLS) {
    method=METHOD_BVLS;
  }
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("ttacfile := %s\n", ttacfile);
    printf("rtacfile := %s\n", rtacfile);
    printf("resfile := %s\n", resfile);
    printf("fitfile := %s\n", fitfile);
    printf("svgfile := %s\n", svgfile);
    printf("method := %s\n", method_str[method]);
    printf("model := %s\n", model_str[model]);
    printf("required_fittime := %g min\n", fitdur);
    printf("weights := %d\n", weights);
    if(!isnan(fixed_k2r)) printf("fixed_k2r := %g\n", fixed_k2r);
    fflush(stdout);
  }
 
 
 
  /*
   *  Read TAC files
   */
  TAC ttac; tacInit(&ttac);
 
  /* Read tissue TAC file */
  if(verbose>1) printf("reading %s\n", ttacfile);
  ret=tacRead(&ttac, ttacfile, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac); return(2);
  }
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(ttac.format));
    printf("tacNr := %d\n", ttac.tacNr);
    printf("sampleNr := %d\n", ttac.sampleNr);
    printf("xunit := %s\n", unitName(ttac.tunit));
    printf("yunit := %s\n", unitName(ttac.cunit));
    printf("midxrange := %g - %g\n", ttac.x[0], ttac.x[ttac.sampleNr-1]);
  }
  if(tacVerifyTimeOrder(&ttac, &status)) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac); return(2);
  }
 
  /* Read reference tissue TAC file */
  if(verbose>1) {printf("reading %s\n", rtacfile); fflush(stdout);}
  int refindex;
  ret=tacReadReference(&ttac, rtacfile, NULL, &refindex, &status);
  if(ret==0) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error)); fflush(stderr);
    tacFree(&ttac); return(2);
  }
  if(ret>1) {
    fprintf(stderr, "Warning: several reference regions found: %s selected.\n",
            ttac.c[refindex].name);
  } else if(verbose>1) {
    printf("reference_region := %s\n", ttac.c[refindex].name); fflush(stderr);
  }
 
  /* Try to fix missing values, if any */
  if(tacFixNaNs(&ttac)) {
    fprintf(stderr, "Error: missing sample(s) in data.\n"); fflush(stderr);
    tacFree(&ttac); return(2);
  }
 
  /* Convert sample times to minutes */
  if(tacXUnitConvert(&ttac, UNIT_MIN, &status)) {
    fprintf(stderr, "Warning: %s\n", errorMsg(status.error)); fflush(stderr);
    // Do not stop in error, assume that user has set units correctly
  }
 
  /* Remove small frame overlaps and gaps */
  if(tacCorrectFrameOverlap(&ttac, &status)) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error)); fflush(stderr);
    tacFree(&ttac); return(2);
  }
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(ttac.format));
    printf("tacNr := %d\n", ttac.tacNr);
    printf("sampleNr := %d\n", ttac.sampleNr);
    printf("xunit := %s\n", unitName(ttac.tunit));
    printf("yunit := %s\n", unitName(ttac.cunit));
    printf("midxrange := %g - %g\n", ttac.x[0], ttac.x[ttac.sampleNr-1]);
  }
 
 
 
  /* Add data weights, if requested */
  if(weights==1) {
    ttac.weighting=WEIGHTING_OFF; 
    for(int i=0; i<ttac.sampleNr; i++) ttac.w[i]=1.0;
  } else if(weights==2) {
    if(tacWByFreq(&ttac, ISOTOPE_UNKNOWN, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&ttac); return(2);
    }
  } else if(!tacIsWeighted(&ttac)) {
    if(verbose>0) fprintf(stderr, "Warning: data is not weighted.\n");
  }
 
  /* Set fit duration */
  int origSampleNr=ttac.sampleNr;
  double starttime=0.0;
  double endtime=fitdur;
  int fitSampleNr=tacFittime(&ttac, &starttime, &endtime, NULL, NULL, &status);
  if(fitSampleNr<0) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac); return(2);
  } else if(fitSampleNr<4) {
    fprintf(stderr, "Error: too few data points for a decent fit.\n");
    tacFree(&ttac); return(2);
  }
  if(verbose>1) {
    printf("starttime := %g\n", starttime);
    printf("endtime := %g\n", endtime);
    //printf("first := %d\n", first);
    //printf("last := %d\n", last);
    printf("fitSampleNr := %d\n", fitSampleNr);
    printf("origSampleNr := %d\n", origSampleNr);
    fflush(stdout);
  }
  fitdur=endtime;
  ttac.sampleNr=fitSampleNr;
 
 
 
  /*
   *  Calculate AUC curves
   */
  TAC auc; tacInit(&auc);
  if(tacInterpolate(&ttac, &ttac, NULL, &auc, NULL, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac); return(3);
  }
  if(verbose>3) {
    printf("\nAUC 0-%g:\n", auc.x[fitSampleNr-1]);
    for(int ti=0; ti<auc.tacNr; ti++)
      printf("  %s\t%g\n", auc.c[ti].name, auc.c[ti].y[fitSampleNr-1]);
    fflush(stdout);
  }
 
  /*
   *  Set model-dependent parameters
   */
  int llsq_n=0;
  switch(model) {
    case MODEL_R1:    llsq_n=2; break;
    case MODEL_R1R2:  llsq_n=3; break;
    default: exit(1);
  }
  if(verbose>2) printf("llsq_n := %d\n", llsq_n);
 
 
  /*
   *  Prepare the room for LLSQ parameters
   */
  if(verbose>1) printf("initializing LLSQ parameter data\n");
  PAR lp; parInit(&lp);
  ret=parAllocateWithTAC(&lp, &ttac, MAX_LLSQ_N, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac);  tacFree(&auc); 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, "datafile", ttacfile, 0, NULL);
    iftPut(&lp.h, "refname", ttac.c[refindex].name, 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(&ttac); tacFree(&auc); 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(&ttac); tacFree(&auc); parFree(&lp); free(llsq_mat); return(5);
      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 (except reference TAC)
     */
    if(verbose>1) printf("solving NNLS\n");
    for(int ti=0; ti<ttac.tacNr; ti++) if(ti!=refindex) {
 
      if(verbose>2 && ttac.tacNr>2) {
        printf("Region %d %s\n", 1+ti, ttac.c[ti].name); fflush(stdout);}
 
      /* Setup data matrix A and vector B */
      for(int mi=0; mi<llsq_m; mi++)
        llsq_b[mi]=ttac.c[ti].y[mi];
      if(model==MODEL_R1) {
        for(int mi=0; mi<llsq_m; mi++) {
          llsq_mat[mi]=ttac.c[refindex].y[mi];                        // RTAC
          llsq_mat[mi+llsq_m]=auc.c[refindex].y[mi]-auc.c[ti].y[mi];  // RTAC - TTAC integral
        }
      } else {
        for(int mi=0; mi<llsq_m; mi++) {
          llsq_mat[mi]=ttac.c[refindex].y[mi];        // RTAC
          llsq_mat[mi+llsq_m]=auc.c[refindex].y[mi];  // RTAC integral
          llsq_mat[mi+2*llsq_m]=-auc.c[ti].y[mi];     // -TTAC integral
        }
      }
      if(verbose>10) {
        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(&ttac)) nnlsWght(llsq_n, llsq_m, llsq_a, llsq_b, ttac.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", ttac.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", ttac.c[ti].name);
      }
      if(verbose>5) {
        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(&ttac);
      lp.r[ti].fitNr=llsq_n;
 
      /* Compute fitted TAC (into auc since it is otherwise not needed) */
      for(int mi=0; mi<llsq_m; mi++) {
        auc.c[ti].y[mi]=0.0;
        for(int ni=0; ni<llsq_n; ni++) auc.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(&ttac);  tacFree(&auc); 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 (except reference TAC)
     */
    if(verbose>1) printf("solving BVLS\n");
    for(int ti=0; ti<ttac.tacNr; ti++) if(ti!=refindex) {
 
      if(verbose>2 && ttac.tacNr>2) {
        printf("Region %d %s\n", 1+ti, ttac.c[ti].name); fflush(stdout);}
 
      /* Setup data matrix A and vector B */
      for(int mi=0; mi<llsq_m; mi++)
        b[mi]=ttac.c[ti].y[mi];
      if(model==MODEL_R1) {
        for(int mi=0; mi<llsq_m; mi++) {
          llsq_mat[mi]=ttac.c[refindex].y[mi];                        // RTAC
          llsq_mat[mi+llsq_m]=auc.c[refindex].y[mi]-auc.c[ti].y[mi];  // RTAC - TTAC integral
        }
      } else {
        for(int mi=0; mi<llsq_m; mi++) {
          llsq_mat[mi]=ttac.c[refindex].y[mi];        // RTAC
          llsq_mat[mi+llsq_m]=auc.c[refindex].y[mi];  // RTAC integral
          llsq_mat[mi+2*llsq_m]=-auc.c[ti].y[mi];     // -TTAC integral
        }
      }
      if(verbose>10) {
        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(&ttac)) llsqWght(llsq_n, llsq_m, NULL, llsq_mat, b, ttac.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;
      if(!isnan(fixed_k2r)) {
        /* Set istate vector to indicate that: */
        istate[llsq_n]=1; // one parameter is bound
        istate[0]=2; // 1+index of that one bound parameter is 2
        istate[1]=1; // 1+index of first active parameter is 1
        istate[2]=3; // 1+index of second active parameter is 3
      }
      /* Set parameter limits */
      if(model==MODEL_R1) {
        bl[0]=0.0; bu[0]=100.0;
        bl[1]=0.0; bu[1]=10.0;
      } else {
        bl[0]=0.0; bu[0]=100.0;
        bl[1]=0.0; bu[1]=10.0;
        bl[2]=0.0; bu[2]=10.0;
        if(!isnan(fixed_k2r)) {
          bl[1]=bu[1]=fixed_k2r;
        }
      }
      /* Set max iterations */
      iterNr=4*llsq_n;
      /* Compute BVLS */
      if(verbose>3) printf("starting BVLS...\n");
      if(isnan(fixed_k2r))
        ret=bvls(0, llsq_m, llsq_n, llsq_mat, b, bl, bu, x, w, act, zz, istate, &iterNr, verbose-3);
      else
        ret=bvls(1, 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", ttac.c[ti].name);
        else fprintf(stderr, "Warning: no BVLS solution for %s\n", ttac.c[ti].name);
        for(int ni=0; ni<llsq_n; ni++) x[ni]=0.0;
        r2=0.0;
      }
      if(verbose>5) {
        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(&ttac);
      lp.r[ti].fitNr=llsq_n;
 
      /* Compute fitted TAC (into auc since it is otherwise not needed) */
      for(int mi=0; mi<llsq_m; mi++) {
        auc.c[ti].y[mi]=0.0;
        for(int ni=0; ni<llsq_n; ni++) auc.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(&ttac); tacFree(&auc); parFree(&lp);
    return(1);
 
  }
 
 
  /*
   *  Print original LLSQ parameters on screen
   */
  if(verbose>2 && lp.tacNr<50) {
    parWrite(&lp, stdout, PAR_FORMAT_TSV_UK, 0, &status);
    fflush(stdout);
  }
 
  /*
   *  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, &ttac, MAX_LLSQ_N+4, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac); tacFree(&auc); parFree(&lp);
    return(21);
  }
  iftFree(&cmpar.h);
  /* Copy titles & file names & parameter values */
  {
    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, "datafile", ttacfile, 0, NULL);
    iftPut(&cmpar.h, "refname", ttac.c[refindex].name, 0, NULL);
    /* Study number */
    if(tacGetHeaderStudynr(&ttac.h, buf, NULL)==TPCERROR_OK)
      iftPut(&cmpar.h, "study_number", buf, 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=lp.r[i].end;
      /* and nr of fitted and saved parameters */
      cmpar.r[i].fitNr=lp.r[i].fitNr;
      /* and WSS */
      cmpar.r[i].wss=lp.r[i].wss;
    }
    /* Set the parameter names and units */
    if(model==MODEL_R1) {
      cmpar.parNr=llsq_n;
      i=0; strcpy(cmpar.n[i].name,"R1"); 
      i++; strcpy(cmpar.n[i].name,"k2"); cmpar.n[i].unit=UNIT_PER_MIN;
    } else {
      cmpar.parNr=llsq_n;
      i=0; strcpy(cmpar.n[i].name,"R1"); 
      i++; strcpy(cmpar.n[i].name,"k2"); cmpar.n[i].unit=UNIT_PER_MIN;
      i++; strcpy(cmpar.n[i].name,"k2r"); cmpar.n[i].unit=UNIT_PER_MIN;
      cmpar.parNr+=2; // two derived parameters
      i++; strcpy(cmpar.n[i].name,"R2"); 
      i++; strcpy(cmpar.n[i].name,"pRatio"); 
    }
    /* Set parameters */
    if(model==MODEL_R1) {
      for(i=0; i<cmpar.tacNr; i++) {
        cmpar.r[i].p[0]=lp.r[i].p[0];
        cmpar.r[i].p[1]=lp.r[i].p[1];
      }
    } else {
      for(i=0; i<cmpar.tacNr; i++) {
        double R1, k2, k2r, R2, pRatio;
        R1=lp.r[i].p[0]; k2=lp.r[i].p[2]; k2r=lp.r[i].p[1]/lp.r[i].p[0];
        R2=k2/k2r; pRatio=R1/R2;
        cmpar.r[i].p[0]=R1;
        cmpar.r[i].p[1]=k2;
        cmpar.r[i].p[2]=k2r;
        cmpar.r[i].p[3]=R2;
        cmpar.r[i].p[4]=pRatio;
      }
    }
  }
  /* Delete reference TAC from parameters */
  parDeleteTAC(&cmpar, refindex);
 
  /*
   *  Print CM parameters on screen
   */
  if(verbose>0 && lp.tacNr<80) parWrite(&cmpar, stdout, PAR_FORMAT_TSV_UK, 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(&ttac); tacFree(&auc); parFree(&lp); parFree(&cmpar); return(22);
    }
    if(verbose>0) printf("Results saved in %s.\n", resfile);
  }
 
  /* Replace 'fitted' (AUC of) reference TAC with original measured TAC */
  for(int i=0; i<auc.sampleNr; i++) auc.c[refindex].y[i]=ttac.c[refindex].y[i];
 
 
  /*
   *  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(&ttac.h, "studynr", 0);
    if(i<0) i=iftFindKey(&ttac.h, "study_number", 0);
    if(i>=0) {strcat(buf, ": "); strcat(buf, ttac.h.item[i].value);}
    ret=tacPlotFitSVG(&ttac, &auc, buf, 0.0, nan(""), 0.0, nan(""), svgfile, &status);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&ttac); tacFree(&auc); 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(&auc, 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(&ttac); tacFree(&auc); parFree(&lp); parFree(&cmpar);
      return(32);
    }
    if(verbose>0) printf("fitted TACs saved in %s.\n", fitfile);
  }
 
 
  tacFree(&ttac); tacFree(&auc); parFree(&lp); parFree(&cmpar);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/