simttm.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 "tpcli.h"
#include "tpccm.h"
#include "tpctacmod.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Simulation of input function time lag using n-compartmental transit-time",
  "model (one-way catenary compartmental model):",
  "  ____        ____        ____              ____       ",
  " |    |  k   |    |  k   |    |  k         |    |  k   ",
  " | Ca | ---> | C1 | ---> | C2 | --->...--> | Cn | ---> ",
  " |____|      |____|      |____|            |____|      ",
  " ",
  " ",
  "Usage: @P [Options] tacfile n k [outputfile] ",
  " ",
  "Options:",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "If file name for output is not given, then the TAC file is modified.",
  "The units of k must correspond the time units of the TAC.",
  "For accurate results, input TAC should have very short sample intervals.",
  "Frame durations are not used, even if available in the TAC file.",
  " ",
  "See also: simdisp, convexpf, fit2dat, simframe, sim_av, tacunit",
  " ",
  "Keywords: TAC, input, simulation, dispersion, time delay",
  0};
/*****************************************************************************/
 
/*****************************************************************************/
/* Turn on the globbing of the command line, since it is disabled by default in
   mingw-w64 (_dowildcard=0); in MinGW32 define _CRT_glob instead, if necessary;
   In Unix&Linux wildcard command line processing is enabled by default. */
/*
#undef _CRT_glob
#define _CRT_glob -1
*/
int _dowildcard = -1;
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int          ai, help=0, version=0, verbose=1;
  char         tacfile[FILENAME_MAX], simfile[FILENAME_MAX];
  double       k=nan("");
  unsigned int n=0;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  tacfile[0]=simfile[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;
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    return(1);
  } else break; // tac name argument may start with '-'
 
  TPCSTATUS status; statusInit(&status);
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  status.verbose=verbose-1;
  
  /* 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) {
    int i;
    if(atoiCheck(argv[ai], &i) || i<1) {
      fprintf(stderr, "Error: invalid n '%s'.\n", argv[ai]); return(1);}
    n=i; ai++;
  }
  if(ai<argc) {
    k=atofVerified(argv[ai]);
    if(!(k>0.0)) {fprintf(stderr, "Error: invalid k '%s'.\n", argv[ai]); return(1);}
    ai++;
  }
  if(ai<argc) {strlcpy(simfile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Is something missing? */
  if(!(k>0.0)) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  if(!simfile[0]) strcpy(simfile, tacfile);
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("tacfile := %s\n", tacfile);
    printf("simfile := %s\n", simfile);
    printf("n := %u\n", n);
    printf("k := %g\n", k);
    fflush(stdout);
  }
 
 
  /*
   *  Read input TAC
   */
  if(verbose>1) {fprintf(stdout, "reading %s\n", tacfile); fflush(stdout);}
  TAC tac; tacInit(&tac);
  if(tacRead(&tac, tacfile, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s (%s)\n", errorMsg(status.error), tacfile);
    tacFree(&tac); return(2);
  }
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(tac.format));
    printf("tacNr := %d\n", tac.tacNr);
    printf("sampleNr := %d\n", tac.sampleNr);
    printf("xunit := %s\n", unitName(tac.tunit));
    printf("yunit := %s\n", unitName(tac.cunit));
    if(tacIsWeighted(&tac)) printf("weighting := yes\n");
    fflush(stdout);
  }
  /* Check for missing sample times */
  if(tacXNaNs(&tac)>0) {
    fprintf(stderr, "Error: missing frame times.\n");
    tacFree(&tac); return(2);
  }
  /* Check for missing concentrations */
  if(tacYNaNs(&tac, -1)>0) {
    fprintf(stderr, "Error: missing concentrations.\n");
    tacFree(&tac); return(2);
  }
  if(tac.sampleNr<3) {
    fprintf(stderr, "Error: too few samples in input data.\n");
    tacFree(&tac); return(2);
  }
  if(tac.tacNr>1 && verbose>0) {
    fprintf(stderr, "Warning: TAC contains %d TACs\n", tac.tacNr);
    fflush(stderr);
  }
  if(tacSortByTime(&tac, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: invalid sample times.\n");
    tacFree(&tac); return(2);
  }
  if(tac.isframe!=0 && verbose>0) {
    if(tacSetX(&tac, &status)!=TPCERROR_OK) { // make sure that frame middle times are set
      fprintf(stderr, "Error: invalid sample times.\n");
      tacFree(&tac); return(2);
    }
    fprintf(stderr, "Warning: frame durations are ignored.\n");
    fflush(stderr);
  }
 
 
  /* 
   *  Simulation
   */
  if(verbose>1) {fprintf(stdout, "simulating\n"); fflush(stdout);}
  for(int ci=0; ci<tac.tacNr; ci++) {
    double sim[tac.sampleNr];
    if(simTTM(tac.x, tac.c[ci].y, tac.sampleNr, k, n, sim)!=0) {
      fprintf(stderr, "Error: invalid data for simulation.\n");
      tacFree(&tac); return(3);
    }
    for(int i=0; i<tac.sampleNr; i++) tac.c[ci].y[i]=sim[i];
  }
 
 
  /*
   *  Write the file
   */
  {
    if(verbose>1) printf("writing simulated data in %s\n", simfile);
    FILE *fp; fp=fopen(simfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing.\n");
      tacFree(&tac); return(11);
    }
    int ret=tacWrite(&tac, fp, TAC_FORMAT_UNKNOWN, 1, &status);
    fclose(fp); tacFree(&tac);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      return(12);
    }
    if(verbose>=0) printf("Simulated data saved in %s\n", simfile);
  }
 
  return(0);
}
/*****************************************************************************/
/*****************************************************************************/