llsqe3.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[] = {
  "Linear fitting of sum of three exponentials function",
  "  f(x) = p1*exp(p2*x) + p3*exp(p4*x) + p5*exp(p6*x) + p7*exp(p8*x)",
  ", where p8=0, to time-activity curves.",
  " ",
  "Usage: @P [options] TAC results",
  " ",
  "Options:",
  " -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 initial linear model are saved in specified file.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: fit_exp, fit_dexp",
  " ",
  "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
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int     ai, help=0, version=0, verbose=1;
  char    tacfile[FILENAME_MAX], resfile[FILENAME_MAX],
          fitfile[FILENAME_MAX], svgfile[FILENAME_MAX], lpfile[FILENAME_MAX];
  int     weights=WEIGHTING_UNKNOWN; // 0=default, 1=no weighting, 2=frequency
  int     ret;
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  tacfile[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=WEIGHTING_OFF; continue;
    } else if(strcasecmp(cptr, "WF")==0) {
      weights=WEIGHTING_ON_FD; 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(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("tacfile := %s\n", tacfile);
    printf("resfile := %s\n", resfile);
    printf("fitfile := %s\n", fitfile);
    printf("svgfile := %s\n", svgfile);
    printf("lpfile := %s\n", lpfile);
    printf("weights := %d\n", weights);
  }
 
 
 
  /*
   *  Read TAC data
   */
  if(verbose>1) printf("reading TAC data\n");
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  TAC tac0; tacInit(&tac0);
  ret=tacRead(&tac0, tacfile, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac0); return(2);
  }
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(tac0.format));
    printf("tacNr := %d\n", tac0.tacNr);
    printf("sampleNr := %d\n", tac0.sampleNr);
    printf("xunit := %s\n", unitName(tac0.tunit));
    printf("yunit := %s\n", unitName(tac0.cunit));
  }
  if(tac0.sampleNr<=MAX_LLSQ_N) {
    fprintf(stderr, "Error: too few samples for fitting.\n");
    tacFree(&tac0); return(2);
  }
 
  /* Add data weights, if requested */
  if(weights!=WEIGHTING_UNKNOWN) tacSetWeights(&tac0, weights, tac0.sampleNr, NULL);
 
  /* Integrate TACs */
  if(verbose>1) printf("integrating TACs\n");
  TAC tac1; tacInit(&tac1); // 1st integral 
  TAC tac2; tacInit(&tac2); // 2nd integral
  TAC tac3; tacInit(&tac3); // 3rd integral
  ret=tacInterpolate(&tac0, &tac0, NULL, &tac1, &tac2, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac0); tacFree(&tac1); tacFree(&tac2); tacFree(&tac3);
  }
  ret=tacDuplicate(&tac2, &tac3);
  if(ret==TPCERROR_OK) {
    for(int i=0; i<tac2.tacNr; i++) {
      ret=liIntegrate(tac2.x, tac2.c[i].y, tac2.sampleNr, tac3.c[i].y, 3, 0);
      if(ret!=TPCERROR_OK) break;
    }
  }
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: cannot make 3rd integral of TAC.\n");
    tacFree(&tac0); tacFree(&tac1); tacFree(&tac2); tacFree(&tac3);
    return(3);
  }
 
 
  /*
   *  Prepare the room for LLSQ parameters
   */
  if(verbose>1) printf("initializing LLSQ parameter data\n");
  PAR lp; parInit(&lp);
  ret=parAllocateWithTAC(&lp, &tac0, MAX_LLSQ_N, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac0); tacFree(&tac1); tacFree(&tac2); tacFree(&tac3);
    return(4);
  }
  iftFree(&lp.h);
  /* Copy titles & file names */
  lp.parNr=MAX_LLSQ_N;
  {
    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, "tacfile", tacfile, 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=tac0.sampleNr;
      /* and nr of fitted parameters */
      lp.r[i].fitNr=MAX_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);
    }
  }
 
 
  /*
   *  Allocate memory required by linear LSQ
   */
  if(verbose>1) printf("allocating memory for LLSQ\n");
  int llsq_n=MAX_LLSQ_N;
  int llsq_m=tac0.sampleNr;
  double *llsq_mat=(double*)malloc((2*llsq_n*llsq_m)*sizeof(double));
  if(llsq_mat==NULL) {
    fprintf(stderr, "Error: cannot allocate memory.\n");
    tacFree(&tac0); tacFree(&tac1); tacFree(&tac2); tacFree(&tac3);
    parFree(&lp);
    return(5);
  }
  double **llsq_a=(double**)malloc(llsq_m*sizeof(double*));
  if(llsq_a==NULL) {
    fprintf(stderr, "Error: cannot allocate memory.\n");
    tacFree(&tac0); tacFree(&tac1); tacFree(&tac2); tacFree(&tac3);
    parFree(&lp);
    return(5);
  }
  for(int mi=0; mi<llsq_m; mi++) llsq_a[mi]=llsq_mat+mi*llsq_n;
  double r2, llsq_b[llsq_m], llsq_r[llsq_n];
  double *matbackup=llsq_mat+llsq_n*llsq_m;
 
  /*
   *  Fit each regional TAC
   */
  int okNr=0;
  for(int ti=0; ti<tac0.tacNr; ti++) {
 
    if(verbose>1 && tac0.tacNr>1) {
      printf("TAC %d %s\n", 1+ti, tac0.c[ti].name); fflush(stdout);}
 
    /* Setup data matrix A and vector B */
    for(int mi=0; mi<llsq_m; mi++)
      llsq_b[mi]=tac0.c[ti].y[mi];
    for(int mi=0; mi<llsq_m; mi++) {
      llsq_a[mi][0]=tac3.c[ti].y[mi];  // TAC 3rd integral
      llsq_a[mi][1]=tac2.c[ti].y[mi];  // TAC 2nd integral
      llsq_a[mi][2]=tac1.c[ti].y[mi];  // TAC 1st integral
      llsq_a[mi][3]=tac0.x[mi]*tac0.x[mi]*tac0.x[mi]; // x^3
      llsq_a[mi][4]=tac0.x[mi]*tac0.x[mi];            // x^2
      llsq_a[mi][5]=tac0.x[mi];                       // x
      llsq_a[mi][6]=1.0;               // Scale 
    }
    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(&tac0)) qrWeight(llsq_n, llsq_m, llsq_a, llsq_b, tac0.w, NULL);
 
    if(verbose>5) {
      printf("\nA matrix and B vector\n");
      for(int mi=0; mi<llsq_m; mi++) {
        for(int ni=0; ni<llsq_n; ni++) printf("\t%g", llsq_a[mi][ni]);
        printf("\t\t%g\n", llsq_b[mi]);
      }
      printf("\n");
    }
 
    /* Compute QR */
    if(verbose>5) printf("starting QR...\n");
    ret=qrLSQ(llsq_a, llsq_b, llsq_r, llsq_m, llsq_n, &r2);
    if(verbose>5) printf("  ... done.\n");
    if(ret==0) {
      for(int ni=0; ni<llsq_n; ni++) lp.r[ti].p[ni]=llsq_r[ni];
      lp.r[ti].wss=r2;
      okNr++;
    } else {
      fprintf(stderr, "Warning: no solution for TAC %d %s\n", 1+ti, tac0.c[ti].name);
      for(int ni=0; ni<llsq_n; ni++) lp.r[ti].p[ni]=0.0;
      lp.r[ti].wss=0.0;
    }
    lp.r[ti].dataNr=tacWSampleNr(&tac0);
    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++) {
      tac1.c[ti].y[mi]=0.0;
      for(int ni=0; ni<llsq_n; ni++) tac1.c[ti].y[mi]+=llsq_r[ni]*llsq_a[mi][ni];
    }
 
  } // next TAC
 
  /* Free allocated memory */
  free(llsq_a); free(llsq_mat);
 
  if(okNr==0) {
    fprintf(stderr, "Error: no solution found for any of the TACs.\n");
    return(10);
  }
 
  /*
   *  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(&tac0); tacFree(&tac1); tacFree(&tac2); tacFree(&tac3);
      parFree(&lp);
      return(11);
    }
    if(verbose>0) printf("Results saved in %s.\n", lpfile);
  }
 
 
  /*
   *  Prepare the room for exponential function parameters
   */
  if(verbose>1) printf("initializing exponential parameter data\n");
  PAR par; parInit(&par);
  ret=parAllocateWithTAC(&par, &tac0, 8, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac0); tacFree(&tac1); tacFree(&tac2); tacFree(&tac3);
    parFree(&lp);
    return(4);
  }
  iftFree(&par.h);
  /* Copy titles & file names */
  par.parNr=8;
  {
    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, "tacfile", tacfile, 0, NULL);
    /* Set current time to results */
    iftPut(&par.h, "analysis_time", ctime_r_int(&t, buf), 0, NULL);
    /* Descriptions for each TAC */
    for(i=0; i<par.tacNr; i++) {
      par.r[i].dataNr=tacWSampleNr(&tac0);
      par.r[i].fitNr=7;
      par.r[i].start=tac0.x[0]; 
      par.r[i].end=tac0.x[tac0.sampleNr-1];
      par.r[i].model=MF_EXP4;
    }
    /* Set the parameter names and units */
    for(i=0; i<par.parNr; i++) {
      sprintf(par.n[i].name, "P%d", 1+i);
    }
    for(int i=0; i<par.parNr; i+=2) {
      par.n[i].unit=tac0.cunit;
      par.n[i+1].unit=unitInverse(tac0.tunit);
    }
  }
 
 
  /*
   *  Solve exponential function parameters from parameters of linear regression
   */
  okNr=0;
  for(int ti=0; ti<tac0.tacNr; ti++) {
 
    if(verbose>1 && tac0.tacNr>1) {
      printf("TAC %d %s\n", 1+ti, tac0.c[ti].name); fflush(stdout);}
 
    for(int i=0; i<par.parNr; i++) par.r[ti].p[i]=0.0;
    par.r[ti].wss=0.0;
 
    /* Calculate the exponential rate constants as the real solutions of cubic equation */ 
    double A=-lp.r[ti].p[2];
    double B=-lp.r[ti].p[1];
    double C=-lp.r[ti].p[0];
    double x1, x2, x3;
    int roots=rootsCubic(A, B, C, &x1, &x2, &x3);
    if(roots==0) {
      fprintf(stderr, "Warning: no solution for TAC %d %s\n", 1+ti, tac0.c[ti].name);
      continue;
    }
    if(verbose>2) {
      printf("  %d roots: %g", roots, x1);
      if(roots>1) printf(" %g", x2);
      if(roots>2) printf(" %g", x3);
      printf("\n");
    }
    /* Write into results */
    par.r[ti].p[3]=x1;
    if(roots>1) par.r[ti].p[5]=x2;
    if(roots>2) par.r[ti].p[7]=x3;
 
    /* Solve the amplitudes with linear regression */
    llsq_m=tac0.sampleNr;
    llsq_n=1+roots;
    llsq_mat=(double*)malloc((llsq_n*llsq_m)*sizeof(double));
    llsq_a=(double**)malloc(llsq_m*sizeof(double*));
    for(int mi=0; mi<llsq_m; mi++) llsq_a[mi]=llsq_mat+mi*llsq_n;
 
    /* Setup data matrix A and vector B */
    /* Must be done separately for each TAC because parameter number may vary */
    for(int mi=0; mi<llsq_m; mi++)
      llsq_b[mi]=tac0.c[ti].y[mi];
    for(int mi=0; mi<llsq_m; mi++)
      llsq_a[mi][0]=1.0;                   // Scale
    for(int mi=0; mi<llsq_m; mi++)
      llsq_a[mi][1]=exp(x1*tac0.x[mi]);    // 1st exponential
    if(llsq_n>2) for(int mi=0; mi<llsq_m; mi++)
      llsq_a[mi][2]=exp(x2*tac0.x[mi]);    // 2nd exponential
    if(llsq_n>3) for(int mi=0; mi<llsq_m; mi++)
        llsq_a[mi][3]=exp(x3*tac0.x[mi]);  // 3rd exponential
 
    /* Apply data weights */
    if(tacIsWeighted(&tac0)) qrWeight(llsq_n, llsq_m, llsq_a, llsq_b, tac0.w, NULL);
 
    if(verbose>5) {
      printf("\nA matrix and B vector\n");
      for(int mi=0; mi<llsq_m; mi++) {
        for(int ni=0; ni<llsq_n; ni++) printf("\t%g", llsq_a[mi][ni]);
        printf("\t\t%g\n", llsq_b[mi]);
      }
      printf("\n");
    }
 
    /* Compute QR */
    if(verbose>5) printf("starting QR...\n");
    ret=qrLSQ(llsq_a, llsq_b, llsq_r, llsq_m, llsq_n, &r2);
    if(verbose>5) printf("  ... done.\n");
    if(ret!=0) {
      fprintf(stderr, "Warning: no solution for TAC %d %s\n", 1+ti, tac0.c[ti].name);
      continue;
    }
    if(verbose>2) {
      printf("  amplitudes: %g %g", llsq_r[0], llsq_r[1]);
      if(llsq_n>2) printf(" %g", llsq_r[2]);
      if(llsq_n>3) printf(" %g", llsq_r[3]);
      printf("\n");
    }
    /* Write into results */
    par.r[ti].p[0]=llsq_r[0];
    par.r[ti].p[2]=llsq_r[1];
    if(llsq_n>2) par.r[ti].p[4]=llsq_r[2];
    if(llsq_n>3) par.r[ti].p[6]=llsq_r[3];
 
    okNr++;
    /* Free allocated memory */
    free(llsq_a); free(llsq_mat);
  }
  if(okNr==0) {
    fprintf(stderr, "Error: no solution found for any of the TACs.\n");
    tacFree(&tac0); tacFree(&tac2); tacFree(&tac2); tacFree(&tac3);
    parFree(&lp); parFree(&par);
    return(10);
  }
 
 
  /*
   *  Print exponential function parameters
   */
  if(verbose>0 && par.tacNr<80) 
    parWrite(&par, stdout, PAR_FORMAT_TSV_UK, 0, &status);
 
  /*
   *  Save exponential function parameters
   */
  {
    if(verbose>1) printf("writing %s\n", resfile);
    par.format=parFormatFromExtension(resfile);
    if(verbose>2) printf("result file format := %s\n", parFormattxt(par.format));
    if(par.format==PAR_FORMAT_UNKNOWN) par.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(&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(&tac0); tacFree(&tac2); tacFree(&tac2); tacFree(&tac3);
      parFree(&lp); parFree(&par);
      return(21);
    }
    if(verbose>0) printf("Results saved in %s.\n", resfile);
  }
 
 
 
  /*
   *  Compute fitted TAC, if necessary
   */
  if(svgfile[0] || fitfile[0]) {
    if(verbose>1) printf("computing fitted TACs\n");
    // into ttac2 since it is otherwise not needed
    for(int ti=0; ti<tac0.tacNr; ti++) {
      for(int mi=0; mi<tac0.sampleNr; mi++)
        tac2.c[ti].y[mi]=par.r[ti].p[0] + par.r[ti].p[2]*exp(par.r[ti].p[3]*tac0.x[mi]) +
                         par.r[ti].p[4]*exp(par.r[ti].p[5]*tac0.x[mi]) + 
                         par.r[ti].p[6]*exp(par.r[ti].p[7]*tac0.x[mi]);
    }
  }
 
 
  /*
   *  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, "3exp");
    i=iftFindKey(&tac0.h, "studynr", 0);
    if(i<0) i=iftFindKey(&tac0.h, "study_number", 0);
    if(i>=0) {strcat(buf, ": "); strcat(buf, tac0.h.item[i].value);}
    ret=tacPlotFitSVG(&tac0, &tac2, buf, 0.0, nan(""), 0.0, nan(""), svgfile, &status);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac0); tacFree(&tac2); tacFree(&tac2); tacFree(&tac3);
      parFree(&lp); parFree(&par);
      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(&tac2, 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(&tac0); tacFree(&tac1); tacFree(&tac2); tacFree(&tac3);
      parFree(&lp);
      return(32);
    }
    if(verbose>0) printf("fitted TACs saved in %s.\n", fitfile);
  }
 
 
  tacFree(&tac0); tacFree(&tac1); tacFree(&tac2); tacFree(&tac3);
  parFree(&lp); parFree(&par);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/