sim_pkcp.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 "tpccm.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Simulation of drug concentration in plasma using pharmacokinetic compartment",
  "models.", 
  "Available models, and model parameters:",
  "  O1CM; one-compartment model with first-order absorption and elimination;",
  "        ED/V1, ka, ke.",
  "  O2CM; two-compartment model with first-order absorption and elimination;",
  "        ED/V1, ka, kd, kr, ke.",
  " ",
  "Usage: @P [options] parfile [simfile]",
  " ",
  "Options:",
  " -model=<O1CM|O2CM>",
  "     Select the model used to simulate the plasma curve; by default the model",
  "     is read from the parameter file.",
/*  " -nr=<nr of samples>",
  "     Number of samples to simulate; by default read from parfile;",
  "     precise simulation requires frequent samples.",
*/
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "To create a template parameter file, do not enter filename for simulated",
  "data, but give the model with option -model.",
  "Simulated time range is determined based on the parameter file.",
  " ",
  "See also: paucinf, sim_3tcm, tacadd, interpol",
  " ",
  "Keywords: plasma, pharmacokinetics",
  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   simfile[FILENAME_MAX], parfile[FILENAME_MAX];
  char  *cptr;
  unsigned int model=0, parNr=0;
  double startTime=nan(""), endTime=nan(""); // hours
  int    ret, pi;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  simfile[0]=parfile[0]=(char)0;
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(!*cptr) continue;
    if(strncasecmp(cptr, "MODEL=", 6)==0 && strlen(cptr)>7) {
      cptr+=6;
      model=modelCodeIndex(cptr); if(model!=0) continue;
      fprintf(stderr, "Error: invalid model selection.\n"); return(1);
    }
    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; use option --help\n");
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    if(model!=0) printf("model := %s\n", modelCode(model));
    printf("parfile := %s\n", parfile);
    printf("simfile := %s\n", simfile);
  }
  
 
 
  /*
   *  Make template parameter file, if no other filenames were given,
   *  and then quit
   */
  if(!simfile[0]) {
    parNr=modelParNr(model);
    if(model==0 || parNr==0) {
      fprintf(stderr, "Error: unknown simulation model; use option -model.\n");
      return(1);
    }
    PAR par; parInit(&par);
    ret=parAllocate(&par, parNr, 3);
    if(ret) {
      fprintf(stderr, "Error: cannot allocate memory for parameters.\n");
      parFree(&par); return(1);
    }
    par.parNr=parNr; par.tacNr=3;
    for(int i=0; i<par.tacNr; i++) {
      sprintf(par.r[i].name, "Cp%d", 1+i);
      par.r[i].model=model;
      par.r[i].start=0.0;
      par.r[i].end=720.0; // 12h
      par.r[i].dataNr=12;
    }
    if(!strcmp(modelCode(model), "O1CM")) {
      strcpy(par.n[0].name, "ED/V1"); par.n[0].unit=UNIT_UMOL_PER_L;
      strcpy(par.n[1].name, "ka"); par.n[1].unit=UNIT_PER_HOUR;
      strcpy(par.n[2].name, "ke"); par.n[2].unit=UNIT_PER_HOUR;
      par.r[0].p[0]=1.0;
      par.r[0].p[1]=0.6;
      par.r[0].p[2]=0.03;
      for(int i=1; i<par.tacNr; i++) {
        par.r[i].p[0]=1.0*par.r[i-1].p[0];
        par.r[i].p[1]=1.0*par.r[i-1].p[1];
        par.r[i].p[2]=2.0*par.r[i-1].p[2];
      }
    } else if(!strcmp(modelCode(model), "O2CM")) {
      strcpy(par.n[0].name, "ED/V1"); par.n[0].unit=UNIT_UMOL_PER_L;
      strcpy(par.n[1].name, "ka"); par.n[1].unit=UNIT_PER_HOUR;
      strcpy(par.n[2].name, "kd"); par.n[2].unit=UNIT_PER_HOUR;
      strcpy(par.n[3].name, "kr"); par.n[3].unit=UNIT_PER_HOUR;
      strcpy(par.n[4].name, "ke"); par.n[4].unit=UNIT_PER_HOUR;
      par.r[0].p[0]=1.0;
      par.r[0].p[1]=0.6;
      par.r[0].p[2]=0.02;
      par.r[0].p[3]=0.02;
      par.r[0].p[4]=0.04;
      for(int i=1; i<par.tacNr; i++) {
        par.r[i].p[0]=1.0*par.r[i-1].p[0];
        par.r[i].p[1]=1.0*par.r[i-1].p[1];
        par.r[i].p[2]=2.0*par.r[i-1].p[2];
        par.r[i].p[3]=2.0*par.r[i-1].p[3];
        par.r[i].p[4]=1.0*par.r[i-1].p[4];
      }
    } else {
      fprintf(stderr, "Error: incompatible simulation model.\n");
      return(1);
    }
    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);
    }
    ret=parWrite(&par, fp, PAR_FORMAT_CSV_UK, 1, &status);
    fclose(fp);
    parFree(&par); return(0);
  }
 
 
  /*
   *  Read parameter file to get parameters for the simulation, and 
   *  set model in case it was given with command-line option.
   */
  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>5) parWrite(&par, stdout, PAR_FORMAT_TSV_UK, 1, NULL);
 
  if(model!=0) for(int i=0; i<par.tacNr; i++) par.r[i].model=model;
  /* Check the validity of model and its parameters */
  for(int i=0; i<par.tacNr; i++) {
    if(verbose>3) 
      printf("model %s for tac %s\n", modelCode(par.r[i].model), par.r[i].name);
    /* check that we got at least as many parameters as the model has */
    if(modelParNr(par.r[i].model)>(unsigned int)par.parNr) {
      fprintf(stderr, "Error: invalid parameters for selected model.\n");
      parFree(&par); return(3);
    }
    /* check that we can find the obligatory model parameters */
    ret=0;
    if(parFindParameter(&par, "ED/V1")<0) ret++;
    if(parFindParameter(&par, "ka")<0) ret++;
    if(!strcmp(modelCode(par.r[i].model), "O1CM")) {
      if(parFindParameter(&par, "ke")<0) ret++;
    } else if(!strcmp(modelCode(par.r[i].model), "O2CM")) {
      if(parFindParameter(&par, "kd")<0) ret++;
      if(parFindParameter(&par, "kr")<0) ret++;
      if(parFindParameter(&par, "ke")<0) ret++;
    } else {
      fprintf(stderr, "Error: invalid model for this simulation.\n");
      if(verbose>1) printf("model := %s\n", modelCode(par.r[i].model));
      parFree(&par); return(3);
    }
    if(ret) {
      fprintf(stderr, "Error: required parameters not available.\n");
      parFree(&par); return(3);
    }
  }
 
#if(0)
  /* If user did not give sampleNr, then try to get that from parameter file */
  if(sampleNr<1) {
    for(int i=0; i<par.tacNr; i++)
      if(par.r[i].dataNr>sampleNr) sampleNr=par.r[i].dataNr;
    /* If not there either, then use 100 as the default */
    if(sampleNr<1) {
      sampleNr=100;
      if(verbose>0) fprintf(stderr, "Warning: unknown sample nr; set to %d.\n", sampleNr);
    }
  }
#endif
  /* Try to get time range from parameter file */
  for(int i=0; i<par.tacNr; i++) {
    if(isnan(startTime) || par.r[i].start<startTime) startTime=par.r[i].start;
    if(isnan(endTime) || par.r[i].end>endTime) endTime=par.r[i].end;
  }
  startTime/=60.0; endTime/=60.0;
  if(isnan(startTime)) startTime=0.0;
  if(isnan(endTime) || endTime<=startTime) endTime=12.0+fabs(startTime);
  if(verbose>1) {
//    printf("sampleNr := %d\n", sampleNr);
    printf("startTime[h] := %g\n", startTime);
    printf("endTime[h] := %g\n", endTime);
  }
 
 
 
  /*
   *  Allocate space for the simulated data;
   *  must start from 0 time, and have frequent sampling.
   */
  /* Try to find the smallest ka */
  double ka=0.0;
  pi=parFindParameter(&par, "ka");
  if(pi>=0) {
    ka=par.r[0].p[pi];
    for(int i=0; i<par.tacNr; i++)
      if(par.r[i].p[pi]<ka) ka=par.r[i].p[pi];
  }
  if(ka<=0.0) ka=1.0;
  /* Set sample nr */
  int isampleNr;
  isampleNr=(int)simSamples(0.0001/ka, 0.01*endTime, endTime, 2, NULL);
  if(verbose>2) {
    printf("isampleNr := %u\n", isampleNr);
  }
  /* Allocate TAC struct for the simulated plasma curve(s) */
  if(verbose>1) 
    printf("allocating memory for %d samples and %d curves\n", 
           isampleNr, par.tacNr);
  TAC sim; tacInit(&sim);
  if(tacAllocateWithPAR(&sim, &par, isampleNr, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    parFree(&par);
    return(4);
  }
/*
  ret=tacAllocate(&sim, isampleNr, par.tacNr);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(ret));
    parFree(&par); tacFree(&sim); return(4);
  }
*/
  sim.sampleNr=isampleNr; 
  sim.tacNr=par.tacNr;
  sim.weighting=WEIGHTING_OFF;
  pi=parFindParameter(&par, "ED/V1"); if(pi>=0) sim.cunit=par.n[pi].unit;
  pi=parFindParameter(&par, "ka");
  if(pi>=0 && par.n[pi].unit==UNIT_PER_HOUR) {
    sim.tunit=UNIT_HOUR;
  }
  sim.format=TAC_FORMAT_PMOD;
  sim.isframe=0;
  /* Set sample times */
  simSamples(0.0001/ka, 0.01*endTime, endTime, 2, sim.x);
  
 
 
  /*
   *  Simulate data for each parameter set
   */
  double gi[isampleNr], bufi[isampleNr];
  for(int i=0; i<par.tacNr; i++) {
    if(verbose>2 && par.tacNr>1) printf("\nSimulation %d\n", 1+i);
    /* Set curve name */
    strcpy(sim.c[i].name, par.r[i].name);
 
    /* Get the parameter values */
    double k1, ka, kd, kr, ke;
 
    k1=parGetParameter(&par, "ED/V1", i);
    ka=parGetParameter(&par, "ka", i);
    ke=parGetParameter(&par, "ke", i);
    kd=parGetParameter(&par, "kd", i); if(isnan(kd)) kd=0.0;
    kr=parGetParameter(&par, "kr", i); if(isnan(kr)) kr=0.0;
    if(isnan(k1) || isnan(ka) || isnan(ke)) {
      fprintf(stderr, "Error: required parameters not available.\n");
      parFree(&par); tacFree(&sim); return(3);
    }
    if(verbose>3) {
      printf("  ED/V1 := %g\n", k1);
      printf("  ka := %g\n", ka);
      printf("  ke := %g\n", ke);
      printf("  kd := %g\n", kd);
      printf("  kr := %g\n", kr);
    }
    /* Integral of input from GI; (1-e^-ka*t)/ka */
    for(int j=0; j<isampleNr; j++)
      gi[j]=(1.0-exp(-ka*sim.x[j]))/ka;
    /* Simulation of drug concentration in plasma; this equals the concentration
       in the central (1st) compartment. */
    double *cp=sim.c[i].y;
    ret=simC2_i(sim.x, gi, isampleNr, k1, ke, kd, kr, bufi, cp, NULL);
    if(ret!=0) {
      fprintf(stderr, "Error: cannot simulate.\n");
      parFree(&par); tacFree(&sim); return(5);
    }
    if(verbose>8)
      for(int j=0; j<isampleNr; j++)
        printf("\t%g\t%g\n", sim.x[j], cp[j]);
  }
  
 
  /*
   *  Save simulated data 
   */
  if(verbose>1) printf("writing %s\n", simfile);
  FILE *fp; fp=fopen(simfile, "w");
  if(fp==NULL) {
    fprintf(stderr, "Error: cannot open file for writing (%s)\n", simfile);
    parFree(&par); tacFree(&sim); return(11);
  }
  ret=tacWrite(&sim, fp, TAC_FORMAT_UNKNOWN, 1, &status);
  fclose(fp); 
  parFree(&par); tacFree(&sim);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error (%d): %s\n", ret, errorMsg(status.error));
    return(12);
  }
  if(verbose>=0) printf("%s saved.\n", simfile);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/