/******************************************************************* pbpk2cpp.cpp Revision 1.1 initial: 02/04/08 ver 1.0 02/16/08 ver 1.1 02/26/08 ver 1.3 todo: put runge back in print report: totalVolume, totalFlow, ... ********************************************************************/ /* Copyright (C) 2008 Doug Neubauer * * This program is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License * as published by the Free Software Foundation. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ #include #include #include #include #include #include #include const char *precisionName = "double"; const char *expectedVersion = "1.3"; const char *eulersName = "eulers"; const char *rungeName = "rungekutta4"; #define INPUTBUFFERSIZE (8192) typedef struct _PlotSubstance{ char *substance ; struct _PlotSubstance *next; int bloodSubstance; } PlotSubstance ; typedef enum { LEVEL0, LEVELGLOBALS, LEVELSUBSTANCE, LEVELCOMPARTMENT, LEVELCOMPARTMENTSUBSTANCE, LEVELDOSAGE, LEVELDOSAGE1, LEVELPLOTTING, } ; typedef enum { // mass balance types UNKNOWNMASSBALANCECOMPARTMENT, GENERICBLOODCOMPARTMENT, GENERICTISSUECOMPARTMENT, LIVERMETABOLISMCOMPARTMENT, } ; typedef enum { COMPARTMENTSUBSTANCENAME, COMPARTMENTSUBSTANCEEQUATION, COMPARTMENTSUBSTANCEMASSBALANCE, } ; typedef enum { PLOTTIMEINDEX, PLOTTITLEINDEX, XAXISINDEX , EXPERIMENTALFILEINDEX , EXPERIMENTALTITLEINDEX , XSCALEINDEX , YAXISINDEX , YMINIMUMINDEX , YMAXIMUMINDEX , YLOGINDEX , XLOGINDEX , PRINTOUTINDEX, NUMBEROFPLOTPARAMETERS } ; typedef struct _{ char *name ; int index; } PlotParameter ; PlotParameter plotParameters[NUMBEROFPLOTPARAMETERS] = { {"plotTime:", PLOTTIMEINDEX}, {"plotTitle:", PLOTTITLEINDEX}, {"xAxis:" , XAXISINDEX }, {"experimentalFile:" , EXPERIMENTALFILEINDEX }, {"experimentalTitle:" , EXPERIMENTALTITLEINDEX }, {"xScale:" , XSCALEINDEX }, {"yAxis:" , YAXISINDEX }, {"yMinimum:" , YMINIMUMINDEX }, {"yMaximum:" , YMAXIMUMINDEX }, {"yLog:" , YLOGINDEX }, {"xLog:" , XLOGINDEX }, {"printout:", PRINTOUTINDEX }, }; char *AddText( char *s ) { char *rt = new char[strlen(s) + 1]; if (rt == 0) exit(1); strcpy(rt,s); return(rt); } class Equation { public: // accessible by all static Equation *CreateNewEquation ( char *text, Equation *pointerHead, char type1 ) { Equation *pointer; Equation *pointer1; pointer = new Equation( text,type1 ); if (pointerHead) { for (pointer1 = pointerHead; pointer1; pointer1 = pointer1->next) { if (pointer1->next == 0) { pointer1->next = pointer; break; } /* if */ } /* forLoop */ return(pointerHead); } /* if */ return(pointer) ; } char *ReadName(){return(name) ;} char *ReadValue(){return(value) ;} char ReadType(){return(type) ;} Equation *ReadNext(){return(next) ;} protected: // accessible by subclass private: // accessible by this class only Equation( char *text, char type1 ) { char buffer[INPUTBUFFERSIZE] ; char *s; char *s1; name = "?" ; value = "?"; type = type1; next = 0 ; strcpy(buffer, text); s1 = strchr(buffer,'=') ; if (s1 == 0) { fprintf(stderr,"bad equation, missing '=': [%s]\n",text); return; } /* if */ s = s1 + 1; s1--; while ((*s1 == ' ') || (*s1 == '\t')) { // scan white reverse *s1 = 0; s1--; if (s1 < buffer) { fprintf(stdout,"bad equation: [%s]\n",text); return ; } /* if */ } /* while */ name = AddText(buffer) ; while ((*s == ' ') || (*s == '\t')) { // scan white s++; } /* while */ s1 = s + strlen(s) - 1; while ((*s1 == ' ') || (*s1 == '\t')) { // scan white reverse *s1 = 0; s1--; if (s1 < buffer) { fprintf(stdout,"bad equation: [%s]\n",text); return ; } /* if */ } /* while */ value = AddText(s) ; } char *name; char *value; char type; Equation *next; }; class Substance { public: // accessible by all Substance(Substance *pointer) { name = "?" ; equationHead = 0; next = pointer ; } void SetName (char *value){ name = value;} void SetEquation (char *text, char type){ equationHead = Equation::CreateNewEquation(text,equationHead,type);} char *ReadName(){return(name) ;} Substance *ReadNext(){return(next) ;} Equation *ReadEquationHead(){return(equationHead) ;} protected: // accessible by subclass private: // accessible by this class only char *name; Equation *equationHead; Substance *next; }; class CompartmentSubstance { public: // accessible by all CompartmentSubstance(CompartmentSubstance *pointer) { name = "?"; equationHead = 0; massBalance = 0; next = pointer; } void AddSubstanceValue(char *text, int command) { switch (command) { case COMPARTMENTSUBSTANCENAME: name = text; break; case COMPARTMENTSUBSTANCEMASSBALANCE: massBalance = text; break; case COMPARTMENTSUBSTANCEEQUATION: equationHead = Equation::CreateNewEquation(text,equationHead,'C'); break; default: break; } /* switch */ } char *ReadName(){return(name) ;} char *ReadMassBalance(){return(massBalance) ;} Equation *ReadEquationHead(){return(equationHead) ;} CompartmentSubstance *ReadNext(){return(next); } protected: // accessible by subclass private: // accessible by this class only char *name; Equation *equationHead; char *massBalance; CompartmentSubstance *next; }; class Compartment { public: // accessible by all Compartment(Compartment *pointer) { name = "?" ; equationHead = 0; compartmentSubstanceHead = 0; next = pointer ; } void SetName (char *value){ name = value;} void SetEquation (char *text, char type){ equationHead = Equation::CreateNewEquation(text,equationHead,type);} void SetNewSubstance () { compartmentSubstanceHead = new CompartmentSubstance( compartmentSubstanceHead); } void AddSubstanceValue(char *text, int command) { compartmentSubstanceHead->AddSubstanceValue(text,command); } char *ReadName(){return(name) ;} Compartment *ReadNext(){return(next) ;} CompartmentSubstance *ReadCompartmentSubstanceHead(){return(compartmentSubstanceHead) ;} Equation *ReadEquationHead(){return(equationHead) ;} protected: // accessible by subclass private: // accessible by this class only char *name; Compartment *next; Equation *equationHead; CompartmentSubstance *compartmentSubstanceHead; }; class Dosage { public: // accessible by all Dosage( Dosage *pointer ) { name = "?" ; substance = "?"; compartment = "?"; amount = "?"; startTime = "?"; type = "?"; next = pointer ; } void SetName (char *value){ name = value;} void SetSubstance (char *value){ substance = value;} void SetCompartment (char *value){ compartment = value;} void SetAmount (char *value){ amount = value;} void SetStartTime (char *value){ startTime = value;} void SetType (char *value){ type = value;} char *ReadSubstance(){return(substance) ;} char *ReadCompartment(){return(compartment) ;} char *ReadAmount(){return(amount) ;} char *ReadStartTime(){return(startTime) ;} Dosage *ReadNext(){return(next) ;} protected: // accessible by subclass private: // accessible by this class only char *name; char *substance; char *compartment; char *amount; char *startTime; char *type; Dosage *next; }; class World { public: // accessible by all World(){InitializeWorld();} void World::ProcessInput (); void World::PrintPlotHeader (char *filename, char *outputFilename); void World::GenerateSimulationCore (); protected: // accessible by subclass private: // accessible by this class only char *ExtractPlotCompartment(char *text) { static char buffer[INPUTBUFFERSIZE]; char *s; strcpy(buffer, text) ; s = strchr(buffer,'.') ; if (s == 0) { fprintf(stderr,"ExtractPlotCompartment: missing '.' [%s]\n",text); exit(1); } /* if */ *s = 0 ; return(buffer) ; } char *ExtractPlotSubstance(char *text) { static char buffer[INPUTBUFFERSIZE]; char *s; strcpy(buffer, text) ; s = strchr(buffer,'.') ; if (s == 0) { fprintf(stderr,"ExtractPlotSubstance: missing '.' [%s]\n",text); exit(1); } /* if */ return(s + 1) ; } void World::InitializeWorld (); char *World::ExtractParameter (char *inputBuffer); char *World::FindParameter (char *inputBuffer, char *parameter); void World::GenerateSimulationGlobals (); void World::GenerateSimulationCompartmentConstants (); char *World::BuildEquation (char *text, Compartment *compartmentPointer, CompartmentSubstance *compartmentSubstancePointer, int command); char *World::BuildVariable (char *variableBuffer, Compartment *compartmentPointerCurrent, CompartmentSubstance *compartmentSubstancePointerCurrent, int command); void World::GenerateSimulationCompartmentSubstanceConstants (); void World::GenerateSimulationSubstanceConstants (); void World::GenerateInjectSubstance (); void World::GeneratePlotPhysiology (); void World::GenerateEulers (); void World::GenerateRungeKutta (); void World::GenerateSubstanceFunctions (); int lineNumber; char *numericalIntegration; Equation *globalsHead; Dosage *dosageHead; Substance *substanceHead; Compartment *compartmentHead; char *plotValues[NUMBEROFPLOTPARAMETERS]; PlotSubstance *plotSubstancePointer; }; // build variable commands (used?) typedef enum { GLOBALCONSTANTSCOMMAND, COMPARTMENTCONSTANTSCOMMAND, COMPARTMENTSUBSTANCESCOMMAND, PLOTPHYSIOLOGYCOMMAND, MASSBALANCECOMMAND, } ; // Equation, variable rules... Warning: Here there be tygers... char *World::BuildVariable (char *variableBuffer, Compartment *compartmentPointerCurrent, CompartmentSubstance *compartmentSubstancePointerCurrent, int command) { static char buffer[INPUTBUFFERSIZE]; char *s; char *s1; char *compartmentName ; char *substanceName; Equation *equationPointer; Compartment *compartmentPointer; CompartmentSubstance *compartmentSubstancePointer; Substance *substancePointer; compartmentName = "?"; substanceName = "?"; if (compartmentPointerCurrent) { compartmentName = compartmentPointerCurrent->ReadName(); } /* if */ if (compartmentSubstancePointerCurrent) { substanceName = compartmentSubstancePointerCurrent->ReadName(); } /* if */ if (stricmp(variableBuffer, "amount") == 0) { sprintf(buffer,"%s_%s_amount",compartmentName,substanceName); return(buffer) ; } /* if */ s = strchr(variableBuffer,'.') ; if (s == 0) { // variable: step1: search local variables, step2: search compartment variables, step3: search substance variables (todo) , step4: error if (compartmentSubstancePointerCurrent) { for (equationPointer = compartmentSubstancePointerCurrent->ReadEquationHead(); equationPointer; equationPointer = equationPointer->ReadNext()) { if (stricmp(variableBuffer, equationPointer->ReadName()) == 0) { sprintf(buffer,"%s_%s_%s",compartmentName,substanceName,equationPointer->ReadName()); return(buffer) ; } /* if */ } /* forLoop */ } /* if */ if (compartmentPointerCurrent) { for (equationPointer = compartmentPointerCurrent->ReadEquationHead(); equationPointer; equationPointer = equationPointer->ReadNext()) { if (stricmp(variableBuffer, equationPointer->ReadName()) == 0) { sprintf(buffer,"%s_%s",compartmentName,equationPointer->ReadName()); return(buffer) ; } /* if */ } /* forLoop */ } /* if */ for (substancePointer = substanceHead; substancePointer; substancePointer = substancePointer->ReadNext()) { if (stricmp(substancePointer->ReadName(), substanceName) == 0) { for (equationPointer = substancePointer->ReadEquationHead(); equationPointer; equationPointer = equationPointer->ReadNext()) { if (stricmp(variableBuffer, equationPointer->ReadName()) == 0) { sprintf(buffer,"%s_%s",substanceName,equationPointer->ReadName()); return(buffer) ; } /* if */ } /* forLoop */ break; } /* if */ } /* forLoop */ fprintf(stderr,"couldn't find variable: [%s]\n",variableBuffer); return("????") ; } /* if */ s1 = strchr(s + 1,'.') ; if (s1) { // variable.variable.variable ... this ones simple *s = '_' ; *s1 = '_' ; strcpy(buffer, variableBuffer); return(buffer) ; } /* if */ *s = 0; s++; if (stricmp(variableBuffer, "global") == 0) { // global.variable (special case) sprintf(buffer,"%s_%s",variableBuffer,s); return(buffer) ; } /* if */ // printf("test.24 var:[%s] s:[%s] compartmentPointer:%x compartmentSubstancePointer:%x\n",variableBuffer,s,compartmentPointer,compartmentSubstancePointer); // variable.variable if (compartmentPointerCurrent) { for (compartmentSubstancePointer = compartmentPointerCurrent->ReadCompartmentSubstanceHead(); compartmentSubstancePointer; compartmentSubstancePointer = compartmentSubstancePointer->ReadNext()) { if (stricmp(variableBuffer, compartmentSubstancePointer->ReadName()) == 0) { sprintf(buffer,"%s_%s_%s",compartmentName,compartmentSubstancePointer->ReadName(), s); return(buffer) ; } /* if */ } /* forLoop */ } /* if */ for (compartmentPointer = compartmentHead; compartmentPointer; compartmentPointer = compartmentPointer->ReadNext()) { if (stricmp(variableBuffer, compartmentPointer->ReadName()) == 0) { for (compartmentSubstancePointer = compartmentPointer->ReadCompartmentSubstanceHead(); compartmentSubstancePointer; compartmentSubstancePointer = compartmentSubstancePointer->ReadNext()) { if (stricmp(s, compartmentSubstancePointer->ReadName()) == 0) { // compartment.substance (special case) sprintf(buffer,"%s_%s_amount",compartmentPointer->ReadName(),compartmentSubstancePointer->ReadName()); return(buffer) ; } /* if */ } /* forLoop */ sprintf(buffer,"%s_%s_%s",compartmentPointer->ReadName(),substanceName, s); return(buffer) ; } /* if */ } /* forLoop */ fprintf(stderr,"couldn't find variable.variable: [%s]\n",variableBuffer); return("????") ; } char *World::BuildEquation (char *text, Compartment *compartmentPointer, CompartmentSubstance *compartmentSubstancePointer, int command) { static char buffer[INPUTBUFFERSIZE]; char variableBuffer[INPUTBUFFERSIZE]; char *s; int i,j,k; char ch; for (i = 0, j = 0; i < INPUTBUFFERSIZE; ) { ch = text[j] ; if (ch == 0) { break; } /* if */ if (isalpha(ch)) { if (ch == 'e') { if (j > 0) { if ((text[j - 1] >= '0') && (text[j - 1] <= '9')) { buffer[i++] = ch ; // special case: scientific notation: ie: 1.2e-4 j++; continue; } /* if */ } /* if */ } /* if */ for (k = 0; k < INPUTBUFFERSIZE; ) { ch = text[j] ; if ((isalnum(ch)) || (ch == '.')) { // printf("1 ch:[%c]\n",ch); variableBuffer[k++] = ch; j++; } /* if */ else { break; } /* else */ } /* kForLoop */ variableBuffer[k] = 0; // printf("test.1 variable[%s] text[%s]\n",variableBuffer,text); s = BuildVariable(variableBuffer, compartmentPointer, compartmentSubstancePointer,command); // printf("test.2 result: [%s]\n",s); buffer[i] = 0 ; strcat(buffer,s) ; i = strlen(buffer); continue; } /* if */ buffer[i++] = ch ; j++; } /* iForLoop */ buffer[i] = 0 ; return(buffer) ; } void World::GenerateSimulationGlobals () { Equation *equationPointer ; printf("// global constants\n"); for (equationPointer = globalsHead; equationPointer; equationPointer = equationPointer->ReadNext()) { if (equationPointer->ReadType() == 'C') { printf("const %s global_%s = %s ;\n", precisionName, equationPointer->ReadName(), BuildEquation( equationPointer->ReadValue(), 0, 0, GLOBALCONSTANTSCOMMAND ) ); } /* if */ } /* forLoop */ printf("\n"); } void World::GenerateSimulationSubstanceConstants () { Substance *substancePointer; Equation *equationPointer ; printf("// substance constants\n"); for (substancePointer = substanceHead; substancePointer; substancePointer = substancePointer->ReadNext()) { for (equationPointer = substancePointer->ReadEquationHead(); equationPointer; equationPointer = equationPointer->ReadNext()) { if (equationPointer->ReadType() == 'C') { printf("const %s %s_%s = %s ;\n", precisionName, substancePointer->ReadName(), equationPointer->ReadName(), BuildEquation( equationPointer->ReadValue(), 0, 0, COMPARTMENTCONSTANTSCOMMAND ) ); } /* if */ } /* forLoop */ } /* forLoop */ printf("\n"); } void World::GenerateSimulationCompartmentConstants () { Compartment *compartmentPointer; Equation *equationPointer ; printf("// compartment constants\n"); for (compartmentPointer = compartmentHead; compartmentPointer; compartmentPointer = compartmentPointer->ReadNext()) { for (equationPointer = compartmentPointer->ReadEquationHead(); equationPointer; equationPointer = equationPointer->ReadNext()) { printf("const %s %s_%s = %s ;\n", precisionName, compartmentPointer->ReadName(), equationPointer->ReadName(), BuildEquation( equationPointer->ReadValue(), 0, 0, COMPARTMENTCONSTANTSCOMMAND ) ); } /* forLoop */ } /* forLoop */ printf("\n"); } void World::GenerateSimulationCompartmentSubstanceConstants () { Compartment *compartmentPointer; CompartmentSubstance *compartmentSubstancePointer; Equation *equationPointer ; printf("// compartment substance variables and constants\n"); for (compartmentPointer = compartmentHead; compartmentPointer; compartmentPointer = compartmentPointer->ReadNext()) { for (compartmentSubstancePointer = compartmentPointer->ReadCompartmentSubstanceHead(); compartmentSubstancePointer; compartmentSubstancePointer = compartmentSubstancePointer->ReadNext()) { printf("%s %s_%s_amount = 0.0 ;\n", precisionName, compartmentPointer->ReadName(), compartmentSubstancePointer->ReadName() ); for (equationPointer = compartmentSubstancePointer->ReadEquationHead(); equationPointer; equationPointer = equationPointer->ReadNext()) { printf("const %s %s_%s_%s = %s ;\n", precisionName, compartmentPointer->ReadName(), compartmentSubstancePointer->ReadName(), equationPointer->ReadName(), BuildEquation( equationPointer->ReadValue(), compartmentPointer, compartmentSubstancePointer, COMPARTMENTSUBSTANCESCOMMAND ) ); } /* forLoop */ } /* forLoop */ } /* forLoop */ printf("\n"); } void World::GenerateInjectSubstance () { Dosage *dosagePointer; int first; printf("\n"); printf("void InjectSubstance (%s currentTime)\n",precisionName); printf("{\n"); printf(" static int finished = 0;\n"); printf(" if (finished) {\n"); printf(" return; \n"); printf(" } /* if */ \n"); first = 1; for (dosagePointer = dosageHead; dosagePointer; dosagePointer = dosagePointer->ReadNext()) { printf(" if (currentTime >= %s) {\n", dosagePointer->ReadStartTime()); printf(" %s_%s_amount += %s;\n",dosagePointer->ReadCompartment(),dosagePointer->ReadSubstance(), dosagePointer->ReadAmount()); if (first) { first = 0; printf(" finished = 1; \n"); } /* if */ printf(" } /* if */ \n"); printf("\n"); } /* forLoop */ printf("}\n"); } void World::GeneratePlotPhysiology () { Equation *equationPointer; PlotSubstance *plotPointer ; printf("\n"); printf("void PlotPhysiology (%s currentTime)\n",precisionName); printf("{\n"); printf(" static double nextPlot = 0.0;\n"); printf(" if (currentTime >= nextPlot) {\n" ); printf(" nextPlot += %s; \n",plotValues[PLOTTIMEINDEX]); if (strcmp(plotValues[XSCALEINDEX], "1.0") == 0) { printf(" printf(\"%%g \",currentTime);\n"); } /* if */ else { printf(" printf(\"%%g \",currentTime / %s);\n",plotValues[XSCALEINDEX]); } /* else */ for (plotPointer = plotSubstancePointer; plotPointer; plotPointer = plotPointer->next) { if (plotPointer->bloodSubstance) { printf(" printf(\"%%g \",%s_%s_amount / (%s_volume * %s_blood2PlasmaRatio )); // convert Cbloodtotal to Cplasma\n",ExtractPlotCompartment(plotPointer->substance), ExtractPlotSubstance(plotPointer->substance),ExtractPlotCompartment(plotPointer->substance), ExtractPlotSubstance(plotPointer->substance)); } /* if */ else { printf(" printf(\"%%g \",%s_%s_amount / %s_volume);\n",ExtractPlotCompartment(plotPointer->substance), ExtractPlotSubstance(plotPointer->substance),ExtractPlotCompartment(plotPointer->substance)); } /* else */ } /* forLoop */ printf(" printf(\"\\n\");\n"); if (strcmp(plotValues[PRINTOUTINDEX], "?") != 0) { equationPointer = Equation::CreateNewEquation (plotValues[PRINTOUTINDEX] , 0, 'C') ; // printf("[%s] [%s] [%s]\n",plotValues[PRINTOUTINDEX], equationPointer->ReadName(), equationPointer->ReadValue()); printf(" printf(\"# %s: %%g\\n\",%s);\n",equationPointer->ReadName(), BuildEquation( equationPointer->ReadValue(), 0, 0, PLOTPHYSIOLOGYCOMMAND ) ); } /* if */ printf(" } /* if */ \n"); printf("}\n"); } void World::GenerateSubstanceFunctions () { } void World::GenerateRungeKutta () { } void World::GenerateEulers () { Compartment *compartmentPointer; CompartmentSubstance *compartmentSubstancePointer; char *compartmentName ; char *substanceName; printf("void ProcessPhysiologyTop (%s currentTime) // Eulers\n",precisionName); printf("{\n"); for (compartmentPointer = compartmentHead; compartmentPointer; compartmentPointer = compartmentPointer->ReadNext()) { for (compartmentSubstancePointer = compartmentPointer->ReadCompartmentSubstanceHead(); compartmentSubstancePointer; compartmentSubstancePointer = compartmentSubstancePointer->ReadNext()) { compartmentName = compartmentPointer->ReadName(); substanceName = compartmentSubstancePointer->ReadName(); printf(" %s %s_%s_dAmount = 0.0;\n",precisionName, compartmentName, substanceName); } /* forLoop */ } /* forLoop */ printf("\n"); for (compartmentPointer = compartmentHead; compartmentPointer; compartmentPointer = compartmentPointer->ReadNext()) { for (compartmentSubstancePointer = compartmentPointer->ReadCompartmentSubstanceHead(); compartmentSubstancePointer; compartmentSubstancePointer = compartmentSubstancePointer->ReadNext()) { compartmentName = compartmentPointer->ReadName(); substanceName = compartmentSubstancePointer->ReadName(); printf(" %s_%s_dAmount = %s ;\n", compartmentName,substanceName, BuildEquation( compartmentSubstancePointer->ReadMassBalance(), compartmentPointer, compartmentSubstancePointer, MASSBALANCECOMMAND )); } /* if */ } /* forLoop */ printf("\n"); for (compartmentPointer = compartmentHead; compartmentPointer; compartmentPointer = compartmentPointer->ReadNext()) { for (compartmentSubstancePointer = compartmentPointer->ReadCompartmentSubstanceHead(); compartmentSubstancePointer; compartmentSubstancePointer = compartmentSubstancePointer->ReadNext()) { compartmentName = compartmentPointer->ReadName(); substanceName = compartmentSubstancePointer->ReadName(); printf(" %s_%s_amount += %s_%s_dAmount * global_dTime ;\n", compartmentName, substanceName , compartmentName, substanceName); } /* forLoop */ } /* forLoop */ printf("}\n"); } void World::GenerateSimulationCore () { struct tm *newtime ; time_t aclock; time(&aclock); newtime = localtime(&aclock); printf("\n"); printf("// (WARNING: DO NOT EDIT) Machine generated by pbpk2cpp.cpp on %s\n", asctime(newtime) ); printf("// Title: %s\n",plotValues[PLOTTITLEINDEX]); printf("\n"); printf("#include \n"); printf("#include \n"); printf("#include \n"); printf("#include \n"); printf("#include \n"); printf("\n"); GenerateSimulationGlobals(); GenerateSimulationSubstanceConstants(); GenerateSimulationCompartmentConstants(); GenerateSimulationCompartmentSubstanceConstants(); GenerateInjectSubstance (); GeneratePlotPhysiology (); printf("\n"); if (stricmp(numericalIntegration, rungeName) == 0) { GenerateSubstanceFunctions(); GenerateRungeKutta(); } /* if */ else { GenerateEulers(); } /* else */ printf("\n"); printf("\n"); printf("void RunSimulation ()\n"); printf("{\n"); printf(" %s currentTime;\n",precisionName); printf("\n"); printf(" for (currentTime = global_startTime; currentTime < global_stopTime; currentTime += global_dTime) {\n"); printf(" InjectSubstance(currentTime) ;\n"); printf(" ProcessPhysiologyTop(currentTime) ;\n"); printf(" PlotPhysiology(currentTime) ;\n"); printf(" } /* currentTimeForLoop */ \n"); printf("}\n"); printf("\n"); printf("void main ()\n"); printf("{\n"); printf(" RunSimulation();\n"); printf("}\n"); } void World::InitializeWorld () { int i; numericalIntegration = "?"; lineNumber = 0; globalsHead = 0 ; dosageHead = 0; substanceHead = 0; compartmentHead = 0; for (i = 0; i < NUMBEROFPLOTPARAMETERS; i++) { plotValues[ i ] = "?"; } /* iForLoop */ } void World::PrintPlotHeader (char *filename, char *outputFilename) { static char buffer[INPUTBUFFERSIZE]; FILE *filePointer ; char *s; PlotSubstance *pointer; int j; filePointer = fopen(filename,"w") ; if (filePointer == 0) { fprintf(stderr,"can't open: %s\n",filename); exit(1); } /* if */ fprintf(filePointer,"set autoscale\n"); fprintf(filePointer,"set yrange [%s : %s]\n",plotValues[YMINIMUMINDEX], plotValues[YMAXIMUMINDEX]); fprintf(filePointer,"set grid\n"); fprintf(filePointer,"unset log\n"); fprintf(filePointer,"unset label\n"); fprintf(filePointer,"set xtics auto\n"); fprintf(filePointer,"set ytics auto\n"); fprintf(filePointer,"set mxtics 5\n"); fprintf(filePointer,"set mytics 5\n"); if (stricmp(plotValues[YLOGINDEX], "yes") == 0) { fprintf(filePointer,"set logscale y\n"); } /* if */ if (stricmp(plotValues[XLOGINDEX], "yes") == 0) { fprintf(filePointer,"set logscale x\n"); } /* if */ fprintf(filePointer,"set title \"%s\"\n",plotValues[PLOTTITLEINDEX]); fprintf(filePointer,"set ylabel \"%s\"\n", plotValues[YAXISINDEX]); fprintf(filePointer,"set xlabel \"%s\"\n", plotValues[XAXISINDEX]); fprintf(filePointer,"set terminal png\n"); sprintf(buffer,"%s.png",outputFilename); fprintf(filePointer,"set output \"%s\"\n",buffer); s = ""; sprintf(buffer,"%s.dat",outputFilename); fprintf(filePointer,"plot "); for (j = 0, pointer = plotSubstancePointer; pointer; pointer = pointer->next,j++) { fprintf(filePointer,"%s\"%s\" using 1:%d title '%s' with lines",s,buffer,j + 2,pointer->substance); s = ", "; } /* forLoop */ if (strcmp(plotValues[EXPERIMENTALFILEINDEX] , "?") != 0) { fprintf(filePointer,"%s\"%s\" using 1:2 title '%s' with points 3",s,plotValues[EXPERIMENTALFILEINDEX], plotValues[EXPERIMENTALTITLEINDEX]); } /* if */ fprintf(filePointer,"\n"); fclose(filePointer); } char *World::ExtractParameter (char *inputBuffer) { static char buffer[INPUTBUFFERSIZE]; char *s; char *s1; strcpy(buffer, inputBuffer); s1 = strchr(buffer,';') ; if (s1 == 0) { fprintf(stdout,"missing ';' at line %d\n",lineNumber); return(0) ; } /* if */ s = buffer ; while ((*s == ' ') || (*s == '\t')) { // scan white s++; } /* while */ while ((*s1 == ' ') || (*s1 == '\t') || (*s1 == ';')) { // scan white reverse s1--; if (s1 < s) { fprintf(stdout,"missing parameter at line %d\n",lineNumber); return(0) ; } /* if */ } /* while */ s1++ ; *s1 = 0 ; return(s) ; } char *World::FindParameter (char *inputBuffer, char *parameter) { int length = strlen(parameter) ; char *s; if (strnicmp(inputBuffer, parameter, length) == 0) { s = ExtractParameter(inputBuffer + length) ; return(s) ; } /* if */ return(0) ; } void World::ProcessInput () { char inputBuffer[INPUTBUFFERSIZE]; char *s; char *versionNumber = "?"; int stateValue = LEVEL0; int i; PlotSubstance *temporaryPointer; while (1) { if (gets(inputBuffer) == 0){ break; } lineNumber++; if (inputBuffer[0] == 0) { continue; } /* if */ if (inputBuffer[0] == ' ') { continue; } /* if */ if ((inputBuffer[0] == '/') && (inputBuffer[1] == '/') ) { continue; } /* if */ switch (stateValue) { case LEVEL0: if ((s = FindParameter(inputBuffer, "Pbpk2CppVersion:")) != 0) { versionNumber = AddText( s ) ; break; } /* if */ if (strnicmp(inputBuffer, "Global {", 11) == 0) { stateValue = LEVELGLOBALS; break; } /* if */ if (strnicmp(inputBuffer, "Compartment {", 13) == 0) { stateValue = LEVELCOMPARTMENT; compartmentHead = new Compartment(compartmentHead); break; } /* if */ if (strnicmp(inputBuffer, "Substance {", 11) == 0) { stateValue = LEVELSUBSTANCE; substanceHead = new Substance(substanceHead); break; } /* if */ if (strnicmp(inputBuffer, "Dosage {", 11) == 0) { stateValue = LEVELDOSAGE; break; } /* if */ if (strnicmp(inputBuffer, "Plotting {", 10) == 0) { stateValue = LEVELPLOTTING; break; } /* if */ fprintf(stdout,"unknown command: [%s] at line %d \n",inputBuffer,lineNumber); break; case LEVELCOMPARTMENTSUBSTANCE: if (inputBuffer[0] == ']') { stateValue = LEVELCOMPARTMENT; break; } /* if */ if ((s = FindParameter(inputBuffer, "substance:")) != 0) { compartmentHead->AddSubstanceValue(AddText( s ), COMPARTMENTSUBSTANCENAME) ; break; } /* if */ if ((s = FindParameter(inputBuffer, "constant:")) != 0) { compartmentHead->AddSubstanceValue ( AddText(s) , COMPARTMENTSUBSTANCEEQUATION) ; break; } /* if */ if ((s = FindParameter(inputBuffer, "massBalance:")) != 0) { compartmentHead->AddSubstanceValue ( AddText(s) , COMPARTMENTSUBSTANCEMASSBALANCE) ; break; } /* if */ fprintf(stdout,"unknown compartment substance parameter: [%s] at line %d \n",inputBuffer,lineNumber); break; case LEVELCOMPARTMENT: if (inputBuffer[0] == '}') { stateValue = LEVEL0; break; } /* if */ if (inputBuffer[0] == '[') { stateValue = LEVELCOMPARTMENTSUBSTANCE; compartmentHead->SetNewSubstance(); break; } /* if */ if ((s = FindParameter(inputBuffer, "name:")) != 0) { compartmentHead->SetName( AddText( s ) ); break; } /* if */ if ((s = FindParameter(inputBuffer, "constant:")) != 0) { compartmentHead->SetEquation ( AddText( s ) , 'C') ; break; } /* if */ fprintf(stdout,"unknown compartment parameter: [%s] at line %d \n",inputBuffer,lineNumber); break; case LEVELSUBSTANCE: if (inputBuffer[0] == '}') { stateValue = LEVEL0; break; } /* if */ if ((s = FindParameter(inputBuffer, "name:")) != 0) { substanceHead->SetName( AddText( s ) ); break; } /* if */ if ((s = FindParameter(inputBuffer, "constant:")) != 0) { substanceHead->SetEquation (AddText( s ) , 'C') ; break; } /* if */ if ((s = FindParameter(inputBuffer, "string:")) != 0) { substanceHead->SetEquation (AddText( s ) , 'S') ; break; } /* if */ fprintf(stdout,"unknown substance parameter: [%s] at line %d \n",inputBuffer,lineNumber); break; case LEVELGLOBALS: if (inputBuffer[0] == '}') { stateValue = LEVEL0; break; } /* if */ if ((s = FindParameter(inputBuffer, "constant:")) != 0) { globalsHead = Equation::CreateNewEquation (AddText( s ) , globalsHead, 'C') ; break; } /* if */ if ((s = FindParameter(inputBuffer, "string:")) != 0) { globalsHead = Equation::CreateNewEquation ( AddText( s ) , globalsHead, 'S') ; break; } /* if */ if ((s = FindParameter(inputBuffer, "numericalIntegration:")) != 0) { numericalIntegration = AddText( s ) ; break; } /* if */ fprintf(stdout,"unknown global: [%s] at line %d \n",inputBuffer,lineNumber); break; case LEVELDOSAGE: if (inputBuffer[0] == '}') { stateValue = LEVEL0; break; } /* if */ if (inputBuffer[0] == '[') { stateValue = LEVELDOSAGE1; dosageHead = new Dosage(dosageHead); break; } /* if */ fprintf(stdout,"unknown dosage parameter: [%s] at line %d \n",inputBuffer,lineNumber); break; case LEVELDOSAGE1: if (inputBuffer[0] == ']') { stateValue = LEVELDOSAGE; break; } /* if */ if ((s = FindParameter(inputBuffer, "name:")) != 0) { dosageHead->SetName( AddText( s ) ); break; } /* if */ if ((s = FindParameter(inputBuffer, "substance:")) != 0) { dosageHead->SetSubstance( AddText( s ) ); break; } /* if */ if ((s = FindParameter(inputBuffer, "compartment:")) != 0) { dosageHead->SetCompartment( AddText( s ) ); break; } /* if */ if ((s = FindParameter(inputBuffer, "type:")) != 0) { dosageHead->SetType( AddText( s ) ); break; } /* if */ if ((s = FindParameter(inputBuffer, "amount:")) != 0) { dosageHead->SetAmount( AddText( s ) ); break; } /* if */ if ((s = FindParameter(inputBuffer, "startTime:")) != 0) { dosageHead->SetStartTime( AddText( s ) ); break; } /* if */ fprintf(stdout,"unknown dosage parameter: [%s] at line %d \n",inputBuffer,lineNumber); break; case LEVELPLOTTING: if (inputBuffer[0] == '}') { stateValue = LEVEL0; break; } /* if */ if ((s = FindParameter(inputBuffer, "substance:")) != 0) { temporaryPointer = new PlotSubstance() ; temporaryPointer->substance = AddText( s ) ; temporaryPointer->next = plotSubstancePointer; temporaryPointer->bloodSubstance = 0; plotSubstancePointer = temporaryPointer; break; } /* if */ if ((s = FindParameter(inputBuffer, "substanceBlood:")) != 0) { temporaryPointer = new PlotSubstance() ; temporaryPointer->substance = AddText( s ) ; temporaryPointer->next = plotSubstancePointer; temporaryPointer->bloodSubstance = 1; plotSubstancePointer = temporaryPointer; break; } /* if */ for (i = 0; i < NUMBEROFPLOTPARAMETERS; i++) { if ((s = FindParameter(inputBuffer, plotParameters[i].name )) != 0) { plotValues[ plotParameters[i].index ] = AddText( s ) ; break; } /* if */ } /* iForLoop */ if (i == NUMBEROFPLOTPARAMETERS) { fprintf(stdout,"unknown plotting parameter: [%s] at line %d \n",inputBuffer,lineNumber); } /* if */ break; default: break; } /* switch */ } /* while */ if (strcmp(versionNumber, expectedVersion) != 0) { fprintf(stderr,"FATAL: bad versionNumber: expected:[%s] found:[%s]\n",expectedVersion,versionNumber); exit(1); } /* if */ if ( (stricmp(numericalIntegration, eulersName) != 0) && (stricmp(numericalIntegration, rungeName) != 0) ) { fprintf(stderr,"FATAL: bad numericalIntegration: found:[%s]\n",numericalIntegration); exit(1); } /* if */ } int main() { World process ; process.ProcessInput() ; process.PrintPlotHeader("data.plt", "calc2"); process.GenerateSimulationCore(); return(0) ; }