tacsa.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 "tpcbfm.h"
#include "tpctacmod.h"
#include "tpclinopt.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Compartmental model based spectral analysis of regional PET data.",
  "Basis functions (from 1 to n) are formed with a range of k2 from zero to",
  "a given maximum, with the provided PTAC as input function.",
  "Model Croi(t) = Vb*Cb(t) + Ct[1](t) + Ct[2](t) + ... + Ct[n](t)",
  "is fitted to the provided TTACs.", 
  " ",
  "  ____             _________  ",
  " |    |--K1[i=1]->| Ct[i=1] | ",
  " |    |<-k2[i=1]--|_________| ",
  " | Cp |    .......        ",
  " |    |            _________  ",
  " |    |--K1[i=n]->| Ct[i=n] | ",
  " |    |<-k2[i=n]--|_________| ",
  " |____|                       ",
  " ",
  "Croi = Vb*Cb + K1[i=1]*Ct[i=1] + ... K1[i=n]*Ct[i=n]",
  " ",
  "Usage: @P [options] ptacfile btacfile ttacfile maxk2 parfile [safile]",
  " ",
  "Options:",
  " -n=<Number of basis functions>",
  "     Set the number of basis functions; by default 100, minimum 10.",
  " -mink2=<Value>",
  "     Set minimum value for k2 in units 1/min; by default 0.005 min-1.",
  "     Do not set to zero or very low value, but if you want to add basis",
  "     function with k2=0 (trapping compartment), enter minimum k2 as",
  "     negative value.", 
  " -bf=<filename>",
  "     Basis function curves are written in specified file.",
  " -fit=<Filename>",
  "     Fitted regional TTACs are written in specified file.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Parameters Vb and K1 will be written in parameter file as the weights",
  "of the basis functions (Vb as the weight of BTAC, and K1 as the sum of",
  "weights of other basis functions), and k2 as the weighted average of",
  "k2-derived basis functions.",
  "In the optional SA file, the weights for each k2-derived basis functions",
  "are listed in TAC format.",
  " ",
  "See also: fitdelay, bfmh2o, lhsol, hist2svg",
  " ",
  "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;
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int     ai, help=0, version=0, verbose=1;
  char    ptacfile[FILENAME_MAX], btacfile[FILENAME_MAX], ttacfile[FILENAME_MAX],
          bffile[FILENAME_MAX], parfile[FILENAME_MAX], safile[FILENAME_MAX], fitfile[FILENAME_MAX];
  double  k2max=nan("");
  double  k2min=0.005;
  int     bfNr=100;
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  ptacfile[0]=btacfile[0]=ttacfile[0]=bffile[0]=parfile[0]=safile[0]=fitfile[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, "BF=", 3)==0) {
      strlcpy(bffile, cptr+3, FILENAME_MAX); if(strlen(bffile)>0) continue;
    } else if(strncasecmp(cptr, "FIT=", 4)==0) {
      strlcpy(fitfile, cptr+4, FILENAME_MAX); if(strlen(fitfile)>0) continue;
    } else if(strncasecmp(cptr, "N=", 2)==0) {
      if(atoiCheck(cptr+2, &bfNr)==0 && bfNr>=10) continue;
    } else if(strncasecmp(cptr, "NR=", 3)==0) {
      if(atoiCheck(cptr+3, &bfNr)==0 && bfNr>=10) continue;
    } else if(strncasecmp(cptr, "MINK2=", 6)==0) {
      if(atofCheck(cptr+6, &k2min)==0 && fabs(k2min)>=1.0E-04) 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(ptacfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(btacfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(ttacfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {
    if(atofCheck(argv[ai], &k2max) || !(k2max>1.0E-03)) {
      fprintf(stderr, "Error: invalid k2 maximum '%s'.\n", argv[ai]); return(1);}
    ai++;
  }
  if(ai<argc) strlcpy(parfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(safile, 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(!parfile[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("btacfile := %s\n", btacfile);
    printf("ttacfile := %s\n", ttacfile);
    printf("parfile := %s\n", parfile);
    if(safile[0]) printf("safile := %s\n", safile);
    if(bffile[0]) printf("bffile := %s\n", bffile);
    if(fitfile[0]) printf("fitfile := %s\n", fitfile);
    printf("k2max := %g\n", k2max);
    printf("k2min := %g\n", k2min);
    printf("n := %d\n", bfNr);
  }
 
 
  /*
   *  Read tissue and input data
   */
  if(verbose>1) printf("reading tissue and input data\n");
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  double fitdur=1.0E+10;
  TAC ttac, itac; tacInit(&ttac); tacInit(&itac);
  tacReadModelingData(ttacfile, ptacfile, btacfile, NULL, &fitdur, 0, NULL, &ttac, &itac, &status);
  if(status.error!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac); tacFree(&itac); return(2);
  }
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(ttac.format));
    printf("tacNr := %d\n", ttac.tacNr);
    printf("ttac.sampleNr := %d\n", ttac.sampleNr);
    printf("itac.sampleNr := %d\n", itac.sampleNr);
    printf("xunit := %s\n", unitName(ttac.tunit));
    printf("yunit := %s\n", unitName(ttac.cunit));
    printf("fitdur := %g s\n", fitdur);
  }
  if(ttac.sampleNr<5 || itac.sampleNr<5) {
    fprintf(stderr, "Error: too few data samples.\n");
    tacFree(&ttac); tacFree(&itac); return(2);
  }
  /* Set frame mid times */
  if(ttac.isframe) tacSetX(&ttac, NULL);
  if(itac.isframe) tacSetX(&itac, NULL);
  /* Frame times are in minutes, and k2 limits were requested in units 1/min,
     so no need to convert k2min or k2max. */
 
 
  /* Interpolate BTAC to TTAC sample times; needed for Va correction */
  if(verbose>1) printf("interpolating BTAC to TTAC sample times\n");
  TAC tbtac; tacInit(&tbtac);
  tacInterpolate(&itac, &ttac, &tbtac, NULL, NULL, &status);
  if(status.error!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac); tacFree(&itac); return(2);
  }
 
 
  /*
   *  Calculate the basis functions
   */
  if(verbose>1) printf("calculating basis functions\n");
  TAC bf; tacInit(&bf);
  bfm1TCM(&itac, &ttac, bfNr, -k2min, k2max, 1, &bf, &status);
  if(status.error!=TPCERROR_OK) {
    if(verbose>1) fprintf(stderr, "Error: cannot calculate basis functions.\n");
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac); tacFree(&itac); tacFree(&tbtac); return(3);
  }
  if(verbose>2 && bfNr<=20) {
    printf("\nBF k2 values:\n");
    for(int i=0; i<bf.tacNr; i++) printf("\t%d\t%g\n", 1+i, bf.c[i].size);
  }
  /* Save basis functions if required */
  if(bffile[0]) {
    if(verbose>1) printf("writing %s\n", bffile);
    FILE *fp; fp=fopen(bffile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing (%s)\n", bffile);
      tacFree(&ttac); tacFree(&itac); tacFree(&tbtac); tacFree(&bf); 
      return(3);
    }
    tacWrite(&bf, fp, TAC_FORMAT_PMOD, 1, &status);
    fclose(fp);
    if(status.error) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&ttac); tacFree(&itac); tacFree(&tbtac); tacFree(&bf); 
      return(3);
    }
    if(verbose>0) {printf("basis functions saved in %s.\n", bffile); fflush(stdout);}
  }
 
 
  /*
   *  Prepare the room for parameter result
   */
  if(verbose>1) {printf("initializing result data\n"); fflush(stdout);}
  PAR par; parInit(&par);
  parAllocateWithTAC(&par, &ttac, 3, &status);
  if(status.error!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac); tacFree(&itac); tacFree(&tbtac); tacFree(&bf); 
    return(4);
  }
  /* Copy titles & filenames */
  {
    int i;
    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, "plasmafile", ptacfile, 0, NULL);
    iftPut(&par.h, "bloodfile", btacfile, 0, NULL);
    iftPut(&par.h, "datafile", ttacfile, 0, NULL);
    /* Fit method */
    iftPut(&par.h, "fitmethod", "NNLS", 0, NULL);
    /* Set current time to results */
    iftPut(&par.h, "analysis_time", ctime_r_int(&t, buf), 0, NULL);
    /* Set fit times for each TAC */
    for(i=0; i<par.tacNr; i++) {
      par.r[i].dataNr=ttac.sampleNr;
      par.r[i].start=0.0;
      par.r[i].end=fitdur;
      /* and nr of fitted parameters */
      par.r[i].fitNr=3;
    }
    /* Set the parameter names and units */
    i=0; strcpy(par.n[i].name, "K1"); par.n[i].unit=UNIT_ML_PER_ML_MIN;
    i++; strcpy(par.n[i].name, "k2"); par.n[i].unit=UNIT_PER_MIN;
    i++; strcpy(par.n[i].name, "Vb"); par.n[i].unit=UNIT_ML_PER_ML;
  }
 
 
  /*
   *  Allocate memory for regional SA results
   */
  if(verbose>1) {printf("allocating memory for SA results\n"); fflush(stdout);}
  TAC sa; tacInit(&sa);
  {
    if(tacDuplicate(&ttac, &sa)) {
      fprintf(stderr, "Error: cannot allocate memory for SA results.\n");
      tacFree(&ttac); tacFree(&itac); tacFree(&tbtac); tacFree(&bf); parFree(&par);
      return(6);
    }
    if(tacAllocateMoreSamples(&sa, bfNr-sa.sampleNr)) {
      fprintf(stderr, "Error: cannot allocate memory for SA results.\n");
      tacFree(&ttac); tacFree(&itac); tacFree(&tbtac); tacFree(&bf); parFree(&par); tacFree(&sa);
      return(6);
    }
    sa.sampleNr=bfNr; sa.isframe=0; sa.cunit=sa.tunit=UNIT_UNKNOWN;
    for(int i=0; i<bfNr; i++) sa.x[i]=bf.c[i].size;
  }
 
  /*
   *  Allocate memory for fitted TTACs, if requested
   */
  TAC ftac; tacInit(&ftac);
  if(fitfile[0]) {
    if(verbose>1) {printf("allocating memory for fitted TTACs\n"); fflush(stdout);}
    if(tacDuplicate(&ttac, &ftac)) {
      fprintf(stderr, "Error: cannot allocate memory for fitted TTACs.\n");
      tacFree(&ttac); tacFree(&itac); tacFree(&tbtac); tacFree(&bf); parFree(&par); tacFree(&sa);
      return(6);
    }
  }
 
 
  /*
   *  Allocate memory required by NNLS
   */
  if(verbose>1) {printf("allocating memory for NNLS\n"); fflush(stdout);}
  int llsq_m=ttac.sampleNr;
  int llsq_n=1+bf.tacNr; // Including Vb as the first column and parameter
  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(&itac); tacFree(&tbtac); tacFree(&bf); parFree(&par); tacFree(&sa); tacFree(&ftac);
    return(7);
  }
  double **llsq_a=(double**)malloc(llsq_n*sizeof(double*));
  if(llsq_a==NULL) {
    fprintf(stderr, "Error: cannot allocate memory for NNLS.\n");
    fprintf(stderr, "Error: cannot allocate memory for NNLS.\n");
    tacFree(&ttac); tacFree(&itac); tacFree(&tbtac); tacFree(&bf); parFree(&par); tacFree(&sa); tacFree(&ftac);
    free(llsq_mat);
    return(7);
  }
  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<ttac.tacNr; ti++) {
 
    if(verbose>1 && ttac.tacNr>1) {
      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];    // Measure TTAC into B vector
    for(int mi=0; mi<llsq_m; mi++)
      llsq_mat[mi]=tbtac.c[1].y[mi];  // Interpolated BTAC into first A column
    for(int bi=0; bi<bf.tacNr; bi++)
      for(int mi=0; mi<llsq_m; mi++) // Basis functions into next A columns
        llsq_mat[mi+(1+bi)*llsq_m]=bf.c[bi].y[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];
    /* 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", ttac.c[ti].name);
      for(int ni=0; ni<llsq_n; ni++) llsq_x[ni]=0.0;
      r2=0.0;
    } else if(ret>0) {
      fprintf(stderr, "Warning: maximum iterations reached for %s\n", ttac.c[ti].name);
    }
    /* Copy SA weights, excluding weight for BTAC */
    for(int ni=1; ni<llsq_n; ni++) sa.c[ti].y[ni-1]=llsq_x[ni];
    /* Copy results */
    par.r[ti].p[2]=llsq_x[0]; // Vb
    par.r[ti].p[0]=0.0; // K1
    for(int i=0; i<sa.sampleNr; i++) par.r[ti].p[0]+=sa.c[ti].y[i];
    par.r[ti].p[1]=0.0; // k2
    for(int i=0; i<sa.sampleNr; i++) par.r[ti].p[1]+=sa.c[ti].y[i]*bf.c[i].size;
    if(par.r[ti].p[0]>1.0E-06) par.r[ti].p[1]/=par.r[ti].p[0];
    par.r[ti].wss=r2; // r2
    /* Calculate fitted TTAC, if needed */
    if(fitfile[0]) {
      for(int mi=0; mi<llsq_m; mi++) {
        ftac.c[ti].y[mi]=0.0;
        for(int ni=0; ni<llsq_n; ni++) ftac.c[ti].y[mi]+=llsq_x[ni]*matbackup[mi+ni*llsq_m];
      }
    }
  }
  free(llsq_a); free(llsq_mat);
 
  /*
   *  Print results on screen and/or save in file
   */
  par.format=parFormatFromExtension(parfile);
  if(verbose>2) printf("result file format := %s\n", parFormattxt(par.format));
  if(par.format==PAR_FORMAT_UNKNOWN) par.format=PAR_FORMAT_TSV_UK;
  /* On screen */
  if(verbose>0 && par.tacNr<50) parWrite(&par, stdout, PAR_FORMAT_UNKNOWN, 0, &status);
  /* Save results */
  if(parfile[0]) {
    if(verbose>1) printf("writing %s\n", parfile);
    FILE *fp; fp=fopen(parfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing (%s)\n", parfile);
      tacFree(&ttac); tacFree(&itac); tacFree(&tbtac); tacFree(&bf); parFree(&par); tacFree(&sa); tacFree(&ftac);
      return(11);
    }
    int ret=parWrite(&par, fp, PAR_FORMAT_UNKNOWN, 1, &status);
    fclose(fp);
    //parFree(&par);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&ttac); tacFree(&itac); tacFree(&tbtac); tacFree(&bf); parFree(&par); tacFree(&sa); tacFree(&ftac);
      return(12);
    }
    if(verbose>0) printf("Results saved in %s.\n", parfile);
  }
 
 
  /*
   *  Save SA results (weights for each k2)
   */
  if(safile[0]) {
    if(verbose>1) printf("writing %s\n", safile);
    FILE *fp; fp=fopen(safile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing (%s)\n", safile);
      tacFree(&ttac); tacFree(&itac); tacFree(&tbtac); tacFree(&bf); parFree(&par); tacFree(&sa); tacFree(&ftac);
      return(12);
    }
    int ret=tacWrite(&sa, fp, TAC_FORMAT_PMOD, 1, &status);
    fclose(fp);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&ttac); tacFree(&itac); tacFree(&tbtac); tacFree(&bf); parFree(&par); tacFree(&sa); tacFree(&ftac);
      return(12);
    }
    if(verbose>0) printf("SA results written in %s.\n", safile);
  }
 
 
  /*
   *  Save fitted TTACs, if requested
   */
  if(fitfile[0]) {
    if(verbose>1) printf("writing %s\n", fitfile);
    int ret;
    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(&ftac, 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(&itac); tacFree(&tbtac); tacFree(&bf); parFree(&par); tacFree(&sa); tacFree(&ftac);
      return(13);
    }
    if(verbose>0) printf("fitted TACs saved in %s.\n", fitfile);
  }
 
 
  tacFree(&ttac); tacFree(&itac); tacFree(&tbtac); tacFree(&bf); parFree(&par); tacFree(&sa); tacFree(&ftac);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/