siminput.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 "tpcpar.h"
#include "tpctacmod.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Compartmental model simulation of input function.",
  "Not ready for use.",
  " ",
  "Usage: @P [Options] parfile [outputfile] ",
  " ",
  "Options:",
  " -i=<Interval (s)>",
  "     Sample time interval; too long interval as compared to rate constants",
  "     leads to bias; by default 1 s.",
  " -end=<End time (s)",
  "     Duration of simulation; by default 3600 s.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "To create a template parameter file, do not enter names for other files.",
  "If parameter file does not contain units, then per sec and per mL units",
  "are assumed.",
  " ",
  "Example 1: creating template parameter file",
  "     @P input.par",
  "Example 2: simulating input TAC using parameter file",
  "     @P -i=1 -end=3600 input.par simulated.tac",
  " ",
  "See also: sim_pkcp, sim_av, sim_3tcm, simframe",
  " ",
  "Keywords: input, simulation",
  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 icmparcGetPARs(
  ICMPARC *d,
  PAR *par,
  const int verbose
) {
  if(d==NULL || par==NULL) return(1);
  if(verbose>0) printf("%s()\n", __func__);
 
 
 
 
 
  return(1);
}
/*****************************************************************************/
 
/*****************************************************************************/
int simPTAC(
  double *t,
  const int nr,
  double Ti,
  double kav,
  double k11,
  double k12,
  double ku,
  double *cp
) {
  /* Check for data */
  if(nr<2) return 1;
  if(cp==NULL) return 2;
 
  double cv[nr];
 
  double t_last=0.0; if(t[0]<t_last) t_last=t[0]; 
  double ca_last=0.0;
  double cv_last=0.0;
  double ct1_last=0.0;
  double cai=0.0, cai_last=0.0;
  double cvi=0.0, cvi_last=0.0;
  double ct1i=0.0, ct1i_last=0.0;
  double ct1=0.0;
 
  /* Calculate curves */
  for(int i=0; i<nr; i++) {
    /* delta time / 2 */
    double dt2=0.5*(t[i]-t_last);
    if(dt2<0.0) {
      return(5);
    } else if(dt2>0.0) {
      /* Infusion integral; during infusion: time*level, thereafter (infusion time)*level
         and level is set to 1. */
      double ii; if(t[i]<=0.0) ii=0.0; else if(t[i]<Ti) ii=t[i]; else ii=Ti;
      /* Helpers */
      double bca=cai_last+dt2*ca_last;
      double bcv=cvi_last+dt2*cv_last;
      double bct1=ct1i_last+dt2*ct1_last;
      /* Cv */
      cv[i]=(ii-(kav-k11*k12*kav*dt2*dt2/((1.0+(k11+ku)*dt2)*(1.0+k12*dt2)))*bcv 
             + k11*k12*dt2*bca/((1.0+(k11+ku)*dt2)*(1.0+k12*dt2))
             + k12*bct1/(1.0+k12*dt2)
            ) / (1.0 + dt2*(kav - (k11*k12*kav*dt2*dt2)/((1.0+(k11+ku)*dt2)*(1.0+k12*dt2))));
      cvi+=dt2*(cv[i]+cv_last);
      cp[i]=(kav*cvi - (k11+ku)*bca) / (1 + dt2*(k11+ku));
      cai+=dt2*(cp[i]+ca_last);
      ct1=(k11*cai - k12*bct1) / (1 + k12*dt2);
      ct1i+=dt2*(ct1+ct1_last);
    }
 
 
    /* prepare to the next sample */
    t_last=t[i]; 
    ca_last=cp[i];
    cv_last=cv[i];
    ct1_last=ct1;
    cai_last=cai;
    cvi_last=cvi;
    ct1i_last=ct1i;
  }
 
  return 0;
}
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int      ai, help=0, version=0, verbose=1;
  char     simfile[FILENAME_MAX], parfile[FILENAME_MAX];
  double   interval=1.0, duration=3600.0;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  parfile[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;
    if(strncasecmp(cptr, "I=", 2)==0) {
      interval=atofVerified(cptr+2); if(interval>0.0) continue;
    } else if(strncasecmp(cptr, "END=", 4)==0) {
      duration=atofVerified(cptr+4); if(duration>0.0) 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(parfile, argv[ai++], FILENAME_MAX);
  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(!parfile[0]) {
    fprintf(stderr, "Error: missing parameter file name; use option --help\n");
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("parfile := %s\n", parfile);
    if(simfile[0]) printf("simfile := %s\n", simfile);
    printf("duration := %g\n", duration); 
    printf("interval := %g\n", interval);
    fflush(stdout);
  }
 
 
  /*
   *  Make template parameter file, if no other file names were given, and then quit
   */
  if(!simfile[0]) {
    PAR par; parInit(&par);
    if(parAllocate(&par, 20, 1)) {
      fprintf(stderr, "Error: cannot allocate memory for parameters.\n");
      parFree(&par); return(1);
    }
 
#if(0)
    int pi=0;
    strcpy(par.n[pi].name, "Ti"); par.n[pi].unit=UNIT_SEC;
    pi++; strcpy(par.n[pi].name, "Kav"); par.n[pi].unit=UNIT_PER_SEC;
    pi++; strcpy(par.n[pi].name, "K11"); par.n[pi].unit=UNIT_PER_SEC;
    pi++; strcpy(par.n[pi].name, "K12"); par.n[pi].unit=UNIT_PER_SEC;
    pi++; strcpy(par.n[pi].name, "Ku"); par.n[pi].unit=UNIT_PER_SEC;
    par.parNr=pi+1; 
#endif
 
    int ti=0;
    sprintf(par.r[ti].name, "tac%d", 1+ti);
    par.r[ti].model=0;
    par.tacNr=1;
 
    int pi=0;
    /* Administered parent tracer */
    strcpy(par.n[pi].name, "p_Ti"); par.n[pi].unit=UNIT_SEC; par.r[ti].p[pi]=1.0;
    pi++; strcpy(par.n[pi].name, "p_Tdur"); par.n[pi].unit=UNIT_SEC; par.r[ti].p[pi]=15.0;
    pi++; strcpy(par.n[pi].name, "p_Irate"); par.n[pi].unit=UNIT_SEC_KBQ_PER_ML; par.r[ti].p[pi]=1.0E+04;
    pi++; strcpy(par.n[pi].name, "p_VA"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=1.0;
    pi++; strcpy(par.n[pi].name, "p_AT1"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=0.5;
    pi++; strcpy(par.n[pi].name, "p_T1V"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=0.5;
    pi++; strcpy(par.n[pi].name, "p_AT2"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=0.5;
    pi++; strcpy(par.n[pi].name, "p_T2V"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=0.1;
    pi++; strcpy(par.n[pi].name, "p_AU"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=0.01;
 
    /* Metabolite 1 of parent tracer*/
    pi++; strcpy(par.n[pi].name, "m1_VA"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=1.0;
    pi++; strcpy(par.n[pi].name, "m1_AT1"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=0.5;
    pi++; strcpy(par.n[pi].name, "m1_T1V"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=0.5;
    pi++; strcpy(par.n[pi].name, "m1_AT2"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=0.5;
    pi++; strcpy(par.n[pi].name, "m1_T2V"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=0.1;
    pi++; strcpy(par.n[pi].name, "m1_AU"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=0.05;
    pi++; strcpy(par.n[pi].name, "m1_V"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=0.0;
    pi++; strcpy(par.n[pi].name, "m1_T1"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=0.0;
    pi++; strcpy(par.n[pi].name, "m1_T2"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=0.1;
 
    /* Metabolite 2 of parent tracer */
    pi++; strcpy(par.n[pi].name, "m2_VA"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=1.0;
    pi++; strcpy(par.n[pi].name, "m2_AT1"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=0.5;
    pi++; strcpy(par.n[pi].name, "m2_T1V"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=0.5;
    pi++; strcpy(par.n[pi].name, "m2_AT2"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=0.5;
    pi++; strcpy(par.n[pi].name, "m2_T2V"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=0.1;
    pi++; strcpy(par.n[pi].name, "m2_AU"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=0.02;
    pi++; strcpy(par.n[pi].name, "m2_V"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=0.0;
    pi++; strcpy(par.n[pi].name, "m2_T1"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=0.0;
    pi++; strcpy(par.n[pi].name, "m2_T2"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=0.2;
 
    /* Metabolite 1 of metabolite 2 */
    pi++; strcpy(par.n[pi].name, "m2m1_VA"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=1.0;
    pi++; strcpy(par.n[pi].name, "m2m1_AT1"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=0.5;
    pi++; strcpy(par.n[pi].name, "m2m1_T1V"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=0.5;
    pi++; strcpy(par.n[pi].name, "m2m1_AT2"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=0.5;
    pi++; strcpy(par.n[pi].name, "m2m1_T2V"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=0.1;
    pi++; strcpy(par.n[pi].name, "m2m1_AU"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=0.1;
    pi++; strcpy(par.n[pi].name, "m2m1_V"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=0.0;
    pi++; strcpy(par.n[pi].name, "m2m1_T1"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=0.0;
    pi++; strcpy(par.n[pi].name, "m2m1_T2"); par.n[pi].unit=UNIT_PER_SEC; par.r[ti].p[pi]=0.5;
    par.parNr=pi+1; 
 
#if(0)
    int ti=0;
    sprintf(par.r[ti].name, "tac%d", 1+ti);
    par.r[ti].model=0;
    par.tacNr=1;
    pi=0; par.r[ti].p[pi]=15.0;
    pi++; par.r[ti].p[pi]=0.2;
    pi++; par.r[ti].p[pi]=0.05;
    pi++; par.r[ti].p[pi]=0.01;
    pi++; par.r[ti].p[pi]=0.02;
#endif
 
    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);
      parFree(&par); return(11);
    }
    int ret=parWrite(&par, fp, PAR_FORMAT_TSV_UK, 1, &status);
    fclose(fp); parFree(&par);
    if(ret!=TPCERROR_OK) {fprintf(stderr, "Error: cannot write %s\n", parfile); return(11);}
    return(0);
  }
 
 
  /*
   *  Read parameter file to get parameters for the simulation.
   */
  if(verbose>1) fprintf(stdout, "reading %s\n", parfile);
  PAR par; parInit(&par);
  if(parRead(&par, parfile, &status)) {
    fprintf(stderr, "Error: %s (%s)\n", errorMsg(status.error), parfile);
    parFree(&par); return(3);
  }
  if(verbose>10) parWrite(&par, stdout, PAR_FORMAT_TSV_UK, 1, NULL);
  /* Check that there is just one TAC to simulate */
  if(par.tacNr>0) {
    fprintf(stderr, "Error: simulation of only one TAC is currently supported.\n");
    parFree(&par); return(3);
  }
  /* Set simulation data structure */
  ICMPARC sp; icmparcInit(&sp);
  /* Get and check the model parameters for the parent tracer. */
  {
    int ret=0;
    sp.Ti=parGetParameter(&par, "pTi", 0); 
    if(!(sp.Ti>=0.0)) {fprintf(stderr, "Error: invalid/missing Ti for parent.\n"); ret++;}
    else if(parGetParameterUnit(&par, "pTi")==UNIT_MIN) sp.Ti*=60.0;
    sp.Tdur=parGetParameter(&par, "pTdur", 0); 
    if(!(sp.Tdur>0.0)) {fprintf(stderr, "Error: invalid/missing Tdur for parent.\n"); ret++;}
    else if(parGetParameterUnit(&par, "pTdur")==UNIT_MIN) sp.Tdur*=60.0;
    sp.Irate=parGetParameter(&par, "pIrate", 0); 
    if(!(sp.Irate>0.0)) {fprintf(stderr, "Error: invalid/missing Irate for parent.\n"); ret++;}
    else {
      int u=parGetParameterUnit(&par, "pIrate");
      double f=unitConversionFactor(u, UNIT_SEC_KBQ_PER_ML);
      if(f>0.0) sp.Irate*=f;
    }
    sp.k_BV_BA=parGetParameter(&par, "pVA", 0); 
    if(!(sp.k_BV_BA>0.0)) {fprintf(stderr, "Error: invalid/missing pVA for parent.\n"); ret++;}
    else if(parGetParameterUnit(&par, "pVA")==UNIT_PER_MIN) sp.k_BV_BA/=60.0;
 
    
    if(ret) {parFree(&par); icmparcFree(&sp); return(4);}
  }
 
#if(0)
  {
    int ti, pi, punit;
    /* Check parameter Ti */
    pi=parFindParameter(&par, "Ti");
    if(pi<0) {fprintf(stderr, "Error: missing Ti.\n"); parFree(&par); return(3);}
    punit=par.n[pi].unit;
    if(punit==UNIT_MIN) {
      par.n[pi].unit=UNIT_SEC; for(ti=0; ti<par.tacNr; ti++) par.r[ti].p[pi]*=60.0;
    }
    /* Check parameter Kav */
    pi=parFindParameter(&par, "Kav");
    if(pi<0) {fprintf(stderr, "Error: missing Kav.\n"); parFree(&par); return(3);}
    punit=par.n[pi].unit;
    if(punit==UNIT_PER_MIN) {
      par.n[pi].unit=UNIT_PER_SEC; for(ti=0; ti<par.tacNr; ti++) par.r[ti].p[pi]/=60.0;
    }
    /* Check parameter K11 */
    pi=parFindParameter(&par, "K11");
    if(pi<0) {fprintf(stderr, "Error: missing K11.\n"); parFree(&par); return(3);}
    punit=par.n[pi].unit;
    if(punit==UNIT_PER_MIN) {
      par.n[pi].unit=UNIT_PER_SEC; for(ti=0; ti<par.tacNr; ti++) par.r[ti].p[pi]/=60.0;
    }
    /* Check parameter K12 */
    pi=parFindParameter(&par, "K12");
    if(pi<0) {fprintf(stderr, "Error: missing K12.\n"); parFree(&par); return(3);}
    punit=par.n[pi].unit;
    if(punit==UNIT_PER_MIN) {
      par.n[pi].unit=UNIT_PER_SEC; for(ti=0; ti<par.tacNr; ti++) par.r[ti].p[pi]/=60.0;
    }
    /* Check parameter Ku */
    pi=parFindParameter(&par, "Ku");
    if(pi<0) {fprintf(stderr, "Error: missing Ku.\n"); parFree(&par); return(3);}
    punit=par.n[pi].unit;
    if(punit==UNIT_PER_MIN) {
      par.n[pi].unit=UNIT_PER_SEC; for(ti=0; ti<par.tacNr; ti++) par.r[ti].p[pi]/=60.0;
    }
    if(verbose>1) parWrite(&par, stdout, PAR_FORMAT_TSV_UK, 1, NULL);
  }
#endif
 
 
  /*
   *  Allocate memory for simulated curves
   */
  TAC sim; tacInit(&sim);
  int ret=tacAllocate(&sim, 1+ceil(duration/interval), par.tacNr);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: cannot allocate memory for simulated data.\n");
    parFree(&par); return(5);
  }
  /* Set file format */
  sim.format=TAC_FORMAT_PMOD;
  /* Set sample times */
  sim.tunit=UNIT_SEC;
  sim.isframe=1;
  for(int i=0; i<sim._sampleNr; i++) {
    sim.x1[i]=interval*(double)(i);
    sim.x2[i]=interval*(double)(i+1);
    sim.x[i]=0.5*interval*(double)(2*i+1);
    sim.sampleNr++;
    if(sim.x[i]>=duration) break;
  }
  if(verbose>2) printf("sampleNr := %d\n", sim.sampleNr);
  sim.tacNr=par.tacNr;
 
  /*
   *  Simulate
   */
  for(int ti=0; ti<par.tacNr; ti++) {
    double Ti=parGetParameter(&par, "Ti", ti);
    double kav=parGetParameter(&par, "Kav", ti);
    double k11=parGetParameter(&par, "K11", ti);
    double k12=parGetParameter(&par, "K12", ti);
    double ku=parGetParameter(&par, "Ku", ti);
    if(simPTAC(sim.x, sim.sampleNr, Ti, kav, k11, k12, ku, sim.c[ti].y)!=0) {
      fprintf(stderr, "Error: cannot simulate.\n");
      parFree(&par); tacFree(&sim); return(6);
    }
  }
  parFree(&par);
 
 
  /*
   *  Write the TAC 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 (%s)\n", simfile);
    tacFree(&sim); return(11);
  }
  ret=tacWrite(&sim, fp, TAC_FORMAT_UNKNOWN, 1, &status);
  fclose(fp); tacFree(&sim);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    return(12);
  }
  if(verbose>=0) printf("Simulated TAC saved in %s\n", simfile);
 
  return(0);
}
/*****************************************************************************/
/*****************************************************************************/