/* sbml2xpp.c
This owes much to
SBML Development Group <sbml-team@caltech.edu>
who wrote Matlab translators. I couldn't get it to work
on my old version of Matlab, I just wanted to get XPP code
anyway.
The original code contained in TranslateSBML
was initially developed by:
*
* Sarah Keating
* Science and Technology Research Centre
* University of Hertfordshire
* Hatfield, AL10 9AB
* United Kingdom
*
* http://www.sbml.org
* mailto:sbml-team@caltech.edu
*
I used parts of this code verbatim, but once i sort of understood what
was going on, I modified it rather heavily.
I also used sample code from this group
* Ben Bornstein
* The Systems Biology Markup Language Development Group
* ERATO Kitano Symbiotic Systems Project
* Control and Dynamical Systems, MC 107-81
* California Institute of Technology
* Pasadena, CA, 91125, USA
*
* http://www.cds.caltech.edu/erato
* mailto:sbml-team@caltech.edu
Basic idea:
get all the species names & ids
get all the parameter names & ids
if there is an id, use it. If not use the name
get all rules
get all reactions products, reactants, stoichiometry
if the rule sets a parameter, figure it out.
get all functions
get all events
find all named things which are more than 9 characters
assign them unique characters - xxxx.#
so xxxx is first 4 letters, . is '.' and # is a number
reorder longnames from longest to shortest
search and replace every written string!
write ODE file
write out the rules
write parameters which are constant (fixed=1)
write out initial data for species
write all the reaction formula
write odes
cross fingers...
Notes:
I add pow(x,y) = x^y to keep in compliance with math ML
I also add 2 functions gt(x,y)=x-y and lt(x,y)=y-x
so i dont have to parse the event trigger
*/
#include <stdio.h>
#include <stdlib.h>
#include "sbml/SBMLReader.h"
#include "sbml/SBMLTypes.h"
#define MAXLNAM 1024
#define NP 50 /* max parameters per kinetic rule */
#define NPMAX 500 /* max parameters */
#define MAXR 200 /* max reactions any species is in */
#define MAXPEREV 128
typedef struct {
char *src;
char rep[10];
} LONG_NAMES;
LONG_NAMES long_names[1024];
int lnum=0;
typedef struct {
char *f;
char *tc;
char *v;
} RULE;
RULE *rule;
int Nrule=0;
typedef struct {
char *ev;
char *a[MAXPEREV];
int na;
}EVENT;
EVENT *event;
Nevent=0;
typedef struct {
char *name;
double x0;
char *id;
char *tc;
int c;
int bc;
int r[MAXR];
double s[MAXR];
int nrx;
int rule;
} SPECIES;
SPECIES *X_spec;
int N_spec=0;
/* for each reaction, we have a formula
reactants,products, and stoichiometry
parameters are kept in a separate sturcture.
*/
typedef struct {
char *re[NP];
char *pr[NP];
char *formula;
double spr[NP];
double sre[NP];
int npr;
int nre;
} RXN;
typedef struct {
char *name;
char *formula;
int nargs;
char *arg[32];
} FUN_DEF;
FUN_DEF *funs;
int Nfuns=0;
RXN *rxn;
int Nrxn=0;
typedef struct {
char *name;
char *id;
int fixed;
double z;
char *formula;
int unique;
} PAR;
/* dont always know how many */
PAR par[NPMAX];
int Npar=0;
void GetParameter (Model_t *, unsigned int, unsigned int);
void GetReaction (Model_t *, unsigned int, unsigned int);
void GetSpecies (Model_t *, unsigned int, unsigned int);
void GetListRule (Model_t *, unsigned int, unsigned int);
void GetFunctions(Model_t *m);
void GetEvents(Model_t *m);
add_parameter(char *name, char *id,double z,int f);
add_reaction(int i,char *f,int npr,int nre);
add_rule(int i,char *v,char *f, char *tc);
add_reactant(int i,int j,char *name,double s);
add_product(int i,int j,char *name,double s);
add_species(int i,char *name,char *id,double x0,int bc,int c,char *tc);
int main (int argc, char *argv[])
{
SBMLDocument_t *d;
Model_t *m;
unsigned int level,version;
if (argc != 2)
{
printf("\n usage: s2x <filename>\n\n");
return 1;
}
d = readSBML(argv[1]);
m = SBMLDocument_getModel(d);
SBMLDocument_printWarnings(d, stdout);
SBMLDocument_printErrors (d, stdout);
SBMLDocument_printFatals (d, stdout);
level = SBMLDocument_getLevel(d);
version = SBMLDocument_getVersion(d);
GetSpecies (m, level, version);
GetParameter (m, level, version);
GetListRule (m, level, version);
GetReaction (m, level, version);
GetFunctions(m);
GetEvents(m);
dump_species();
dump_parameters();
dump_rules();
dump_reactions();
dump_funs();
dump_events();
write_ode_file(argv[1]);
}
add_reaction(int i,char *f,int npr,int nre)
{
RXN *r;
r=rxn+i;
r->formula=(char *)malloc(strlen(f)+1);
strcpy(r->formula,f);
r->npr=npr;
r->nre=nre;
}
add_rule(int i,char *v,char *f, char *tc)
{
RULE *r;
r=rule+i;
r->f=(char *)malloc(strlen(f)+1);
strcpy(r->f,f);
r->tc=(char *)malloc(strlen(tc)+1);
strcpy(r->tc,tc);
r->v=(char *)malloc(strlen(v)+1);
strcpy(r->v,v);
check_name_len(r->v);
}
add_parameter(char *name, char *id,double z,int f)
{
int i;
if(!is_blank(name)){
for(i=0;i<Npar;i++)
if(strcmp(name,par[i].name)==0){
printf("Hmm par %d(%s) and %d(%s) are the same\n",Npar,name,i,par[i].name);
return;
}
}
if(!is_blank(id)){
for(i=0;i<Npar;i++)
if(strcmp(id,par[i].id)==0){
printf("Hmm par %d(%s) and %d(%s) are the same\n",Npar,id,i,par[i].id);
return;
}
}
if(strlen(name)>0){
par[Npar].name=(char *)malloc(strlen(name)+1);
strcpy(par[Npar].name,name);
check_name_len(name);
}
if(strlen(id)>0){
par[Npar].id=(char *)malloc(strlen(id)+1);
strcpy(par[Npar].id,id);
check_name_len(id);
}
par[Npar].z=z;
par[Npar].fixed=f;
par[Npar].unique=1;
Npar++;
}
void GetEvents(Model_t *m)
{
const Event_t *e;
const EventAssignment_t *ea;
int n,na;
char *ev;
char *a;
int i,j;
const char *variable;
char *formula;
char big[256];
EVENT *x;
n=Model_getNumEvents(m);
event=(EVENT *)malloc(n*sizeof(EVENT));
Nevent=n;
if(n==0)return;
for(i=0;i<n;i++){
x=event+i;
e=Model_getEvent(m, i);
if ( Event_isSetTrigger(e) )
{
ev= SBML_formulaToString( Event_getTrigger(e) );
x->ev=(char *)malloc(strlen(ev)+1);
strcpy(x->ev,ev);
free(ev);
}
na=Event_getNumEventAssignments(e);
x->na=na;
printf("na=%d\n",na);
for(j=0;j<na;j++){
ea=Event_getEventAssignment(e, j);
if ( EventAssignment_isSetMath(ea) ){
variable = EventAssignment_getVariable(ea);
formula = SBML_formulaToString( EventAssignment_getMath(ea) );
sprintf(big,"%s=%s",variable,formula);
x->a[j]=(char *)malloc(strlen(big)+1);
strcpy(x->a[j],big);
free(formula);
}
}
}
}
void GetFunctions(Model_t *m)
{
const ASTNode_t *math;
const FunctionDefinition_t *fd;
char *formula;
char *sa;
char *name;
FUN_DEF *f;
int i,j;
int n,narg;
n=Model_getNumFunctionDefinitions(m);
if(n==0)return;
Nfuns=n;
funs=(FUN_DEF *)malloc(n*sizeof(FUN_DEF));
for(i=0;i<n;i++){
f=funs+i;
fd=Model_getFunctionDefinition(m, i);
if ( FunctionDefinition_isSetMath(fd) )
{
name=FunctionDefinition_getId(fd);
f->name=(char *)malloc(strlen(name));
check_name_len(name);
strcpy(f->name,name);
math=FunctionDefinition_getMath(fd);
narg=ASTNode_getNumChildren(math)-1;
f->nargs=narg;
if (narg > 0)
{
sa=ASTNode_getName( ASTNode_getLeftChild(math));
f->arg[0]=(char *)malloc(strlen(sa)+1);
strcpy(f->arg[0],sa);
for (j = 1; j<narg; ++j)
{
sa=ASTNode_getName( ASTNode_getChild(math, j) );
f->arg[j]=(char *)malloc(strlen(sa)+1);
strcpy(f->arg[j],sa);
}
}
math = ASTNode_getChild(math, ASTNode_getNumChildren(math) - 1);
formula = SBML_formulaToString(math);
f->formula=(char *)malloc(strlen(formula)+1);
strcpy(f->formula,formula);
free(formula);
}
}
}
/* reaction stuff */
void GetReaction(Model_t *m,
unsigned int level,
unsigned int version )
{
int n=Model_getNumReactions(m);
int i;
Reaction_t *r;
char *formula;
KineticLaw_t *kl;
Parameter_t *p;
SpeciesReference_t *s;
int np,j;
int npr,nre;
char *name,*id;
double value,st;
Nrxn=n;
rxn=(RXN *)malloc(n*sizeof(RXN));
for(i=0;i<n;i++){
r=Model_getReaction(m, i);
if (Reaction_isSetKineticLaw(r))
{
kl = Reaction_getKineticLaw(r);
if ( KineticLaw_isSetMath(kl) )
{
npr=Reaction_getNumProducts(r);
nre=Reaction_getNumReactants(r);
np=KineticLaw_getNumParameters(kl);
formula=KineticLaw_getFormula(kl);
add_reaction(i,formula,npr,nre);
for(j=0;j<np;j++){
p=KineticLaw_getParameter(kl,j);
id=Parameter_getId(p);
name=Parameter_getName(p);
value=Parameter_getValue(p);
if(id==NULL)id=" ";
if(name==NULL)name=" ";
add_parameter(name,id,value,1);
}
for(j=0;j<npr;j++){
s=Reaction_getProduct(r,j);
name=SpeciesReference_getSpecies(s);
st=SpeciesReference_getStoichiometry(s);
add_product(i,j,name,st);
}
for(j=0;j<nre;j++){
s=Reaction_getReactant(r,j);
name=SpeciesReference_getSpecies(s);
st=SpeciesReference_getStoichiometry(s);
add_reactant(i,j,name,st);
}
}
}
}
}
/* i=reactin #, j=reactant number or product number
name is variable and s is stoichiometry */
add_reactant(int i,int j,char *name,double s)
{
RXN *r;
r=rxn+i;
r->re[j]=(char *)malloc(strlen(name)+1);
strcpy(r->re[j],name);
r->sre[j]=s;
}
add_product(int i,int j,char *name,double s)
{
RXN *r;
r=rxn+i;
r->pr[j]=(char *)malloc(strlen(name)+1);
strcpy(r->pr[j],name);
r->spr[j]=s;
}
dump_reactions()
{
int i,j; int npr,nre;
RXN *r;
printf("REACTIONS:\n");
for(i=0;i<Nrxn;i++){
r=rxn+i;
printf("rxn %d: %s \n",i,r->formula);
printf("reactants: ");
npr=r->npr;
nre=r->nre;
for(j=0;j<nre;j++)
printf("%s(%g), ",r->re[j],r->sre[j]);
printf("\nproducts: ");
for(j=0;j<npr;j++)
printf("%s(%g), ",r->pr[j],r->spr[j]);
printf("\n");
}
}
dump_events()
{
int i,j,na;
EVENT *ev;
if(Nevent==0)return;
for(i=0;i<Nevent;i++){
ev=event+i;
printf("global 1 %s {",ev->ev);
na=ev->na;
for(j=0;j<na-1;j++)
printf("%s;",ev->a[j]);
printf("%s}\n",ev->a[na-1]);
}
}
dump_funs()
{
int i,j;
FUN_DEF *f;
for(i=0;i<Nfuns;i++){
f=funs+i;
printf("%s(",f->name);
for(j=0;j<f->nargs;j++)
printf("%s,",f->arg[j]);
printf(")=%s\n",f->formula);
}
}
dump_rules()
{
RULE *r;
int i;
if(Nrule>0)
printf("RULES:\n");
for(i=0;i<Nrule;i++){
r=rule+i;
printf("%s=%s\n",r->v,r->f);
}
}
dump_species()
{
SPECIES *x;
int i;
printf("SPECIES: n i t\n");
for(i=0;i<N_spec;i++){
x=X_spec+i;
printf("%s %s %s %g %d %d \n",x->name,x->id,x->tc,x->x0,x->bc,x->c);
}
}
dump_parameters()
{
int i;
printf("PARAMETERS:\n");
for(i=0;i<Npar;i++)
printf("%d %s %s = %g \n",par[i].fixed,par[i].name,par[i].id,par[i].z);
}
add_species(int i,char *name,char *id,double x0,int bc,int c,char *tc)
{
SPECIES *x;
x=X_spec+i;
if(strlen(name)>0){
x->name=(char *)malloc(strlen(name)+1);
strcpy(x->name,name);
check_name_len(name);
}
else
x->name=NULL;
x->x0=x0;
if(strlen(id)>0){
x->id=(char *)malloc(strlen(id)+1);
check_name_len(id);
strcpy(x->id,id);
}
else
x->id=NULL;
if(strlen(tc)>0){
x->tc=(char *)malloc(strlen(tc)+1);
strcpy(x->tc,tc);
}
else
x->tc=NULL;
x->bc=bc;
x->c=c;
x->nrx=0;
x->rule=0;
}
void GetSpecies ( Model_t *pModel,
unsigned int level,
unsigned int version )
{
int n = Model_getNumSpecies(pModel);
char * pacTypecode;
char * pacName;
char * pacId = NULL;
double dInitialAmount;
int nBoundaryCondition;
int nConstant=0;
int i;
Species_t *pSpecies;
X_spec = (SPECIES *)malloc(n * sizeof (SPECIES));
N_spec=n;
for(i=0;i<n;i++){
pSpecies = Model_getSpecies(pModel, i);
pacTypecode = TypecodeToChar(SBase_getTypeCode(pSpecies));
pacName = Species_getName(pSpecies);
dInitialAmount = Species_getInitialAmount(pSpecies);
nBoundaryCondition = Species_getBoundaryCondition(pSpecies);
if(level==2){
pacId=Species_getId(pSpecies);
nConstant = Species_getConstant(pSpecies);
}
if(pacName==NULL)pacName = " ";
if(pacId==NULL)pacId = " ";
if(pacTypecode==NULL)pacTypecode=" ";
add_species(i,pacName,pacId,dInitialAmount,
nBoundaryCondition,nConstant,pacTypecode);
}
}
void
GetParameter ( Model_t *pModel,
unsigned int level,
unsigned int version )
{
int n = Model_getNumParameters(pModel);
char * pacTypecode;
char * pacName;
char * pacId = NULL;
double dValue;
int nConstant=1;
Parameter_t *pParameter;
int i;
for (i = 0; i < n; i++) {
/* determine the values */
pParameter = Model_getParameter(pModel, i);
/* pacTypecode = TypecodeToChar(SBase_getTypeCode(pParameter)); */
pacName = Parameter_getName(pParameter);
dValue = Parameter_getValue(pParameter);
if (level == 2) {
pacId = Parameter_getId(pParameter);
nConstant = Parameter_getConstant(pParameter);
}
if(pacName==NULL)pacName=" ";
if(pacId==NULL)pacId=" ";
add_parameter(pacName,pacId,dValue,nConstant);
}
}
char *
TypecodeToChar (SBMLTypeCode_t typecode)
{
char * pacTypecode;
switch (typecode)
{
case SBML_COMPARTMENT:
pacTypecode = "SBML_COMPARTMENT";
break;
case SBML_EVENT:
pacTypecode = "SBML_EVENT";
break;
case SBML_EVENT_ASSIGNMENT:
pacTypecode = "SBML_EVENT_ASSIGNMENT";
break;
case SBML_FUNCTION_DEFINITION:
pacTypecode = "SBML_FUNCTION_DEFINITION";
break;
case SBML_KINETIC_LAW:
pacTypecode = "SBML_KINETIC_LAW";
break;
case SBML_MODEL:
pacTypecode = "SBML_MODEL";
break;
case SBML_PARAMETER:
pacTypecode = "SBML_PARAMETER";
break;
case SBML_REACTION:
pacTypecode = "SBML_REACTION";
break;
case SBML_SPECIES:
pacTypecode = "SBML_SPECIES";
break;
case SBML_SPECIES_REFERENCE:
pacTypecode = "SBML_SPECIES_REFERENCE";
break;
case SBML_MODIFIER_SPECIES_REFERENCE:
pacTypecode = "SBML_MODIFIER_SPECIES_REFERENCE";
break;
case SBML_UNIT_DEFINITION:
pacTypecode = "SBML_UNIT_DEFINITION";
break;
case SBML_UNIT:
pacTypecode = "SBML_UNIT";
break;
case SBML_ASSIGNMENT_RULE:
pacTypecode = "SBML_ASSIGNMENT_RULE";
break;
case SBML_ALGEBRAIC_RULE:
pacTypecode = "SBML_ALGEBRAIC_RULE";
break;
case SBML_RATE_RULE:
pacTypecode = "SBML_RATE_RULE";
break;
case SBML_SPECIES_CONCENTRATION_RULE:
pacTypecode = "SBML_SPECIES_CONCENTRATION_RULE";
break;
case SBML_COMPARTMENT_VOLUME_RULE:
pacTypecode = "SBML_COMPARTMENT_VOLUME_RULE";
break;
case SBML_PARAMETER_RULE:
pacTypecode = "SBML_PARAMETER_RULE";
break;
default:
pacTypecode = "ERROR";
break;
}
return pacTypecode;
}
void
GetListRule ( Model_t *pModel,
unsigned int unSBMLLevel,
unsigned int unSBMLVersion )
{
int n = Model_getNumRules(pModel);
/* determine the values */
const char * pacTypecode;
const char * pacFormula = NULL;
const char * pacVariable = NULL;
Rule_t *pRule;
int i;
Nrule=n;
rule=(RULE *)malloc(n*sizeof(RULE));
for (i = 0; i < n; i++) {
/* determine the values */
pRule = Model_getRule(pModel, i);
pacTypecode = TypecodeToChar(SBase_getTypeCode(pRule));
if (unSBMLLevel == 1) {
pacFormula = Rule_getFormula(pRule);
}
else if (unSBMLLevel == 2) {
if (Rule_isSetFormula(pRule) == 1){
pacFormula = Rule_getFormula(pRule);
}
else if (Rule_isSetMath(pRule) == 1) {
Rule_setFormulaFromMath(pRule);
pacFormula = Rule_getFormula(pRule);
}
}
/* values for different types of rules */
switch(SBase_getTypeCode(pRule)) {
case SBML_ASSIGNMENT_RULE:
if (AssignmentRule_isSetVariable((AssignmentRule_t *) pRule) == 1) {
pacVariable = AssignmentRule_getVariable((AssignmentRule_t *) pRule);
}
else {
pacVariable = " ";
}
break;
case SBML_ALGEBRAIC_RULE:
pacVariable = " ";
break;
case SBML_RATE_RULE:
if (RateRule_isSetVariable((RateRule_t *) pRule) == 1) {
pacVariable = RateRule_getVariable((RateRule_t *) pRule);
}
else {
pacVariable = " ";
}
break;
case SBML_SPECIES_CONCENTRATION_RULE:
if (SpeciesConcentrationRule_isSetSpecies((SpeciesConcentrationRule_t *) pRule) == 1) {
pacVariable =
SpeciesConcentrationRule_getSpecies((SpeciesConcentrationRule_t *) pRule);
}
else {
pacVariable = " ";
}
break;
case SBML_COMPARTMENT_VOLUME_RULE:
if (CompartmentVolumeRule_isSetCompartment((CompartmentVolumeRule_t *) pRule) == 1) {
pacVariable = CompartmentVolumeRule_getCompartment((CompartmentVolumeRule_t *) pRule);
}
else {
pacVariable = " ";
}
break;
case SBML_PARAMETER_RULE:
if (ParameterRule_isSetName((ParameterRule_t *)pRule) == 1) {
pacVariable = ParameterRule_getName((ParameterRule_t *) pRule);
}
else {
pacVariable = " ";
}
if (ParameterRule_isSetUnits((ParameterRule_t *) pRule) == 1) {
pacVariable = ParameterRule_getUnits((ParameterRule_t *) pRule);
}
else {
pacVariable = " ";
}
break;
default:
pacVariable = " ";
break;
}
if (pacTypecode == NULL) {
pacTypecode = " ";
}
if (pacFormula == NULL) {
pacFormula = " ";
}
add_rule(i,pacVariable,pacFormula,pacTypecode);
}
}
find_parameter(char *s)
{
int i;
for(i=0;i<Npar;i++){
if((strcmp(s,par[i].name)==0)||(strcmp(s,par[i].id)==0))
return i;
}
return -1;
}
find_species(char *s)
{
SPECIES *sp;
int i;
for(i=0;i<N_spec;i++){
sp=X_spec+i;
if((strcmp(s,sp->name)==0)||(strcmp(s,sp->id)==0))
return i;
}
return -1;
}
mark_rule_pars()
{
int i,j;
RULE *r;
for(i=0;i<Nrule;i++){
r=rule+i;
j=find_parameter(r->v);
if(j>-1){
par[j].fixed=-2;
printf("found %s as %d \n",r->v,j);
}
j=find_species(r->v);
if(j>-1){
(X_spec+j)->rule=1;
printf("found %s as %d \n",r->v,j);
}
}
}
is_blank(char *s)
{
int i;
for(i=0;i<strlen(s);i++)
if(s[i]!=' ')return 0;
return 1;
}
/* this searches all reactions and finds the
partici[ating species and marks them
*/
species_participation()
{
int i,j,k,l;
RXN *r;
SPECIES *s;
for(i=0;i<Nrxn;i++){
r=rxn+i;
/* reactant */
for(j=0;j<r->nre;j++){
k=find_species(r->re[j]);
if(k==-1){
printf("WARNING: species %s not found \n",r->re[j]);
continue;
}
s=X_spec+k;
l=s->nrx;
s->r[l]=i;
s->s[l]=-r->sre[j]; /* change sign for reactant */
l++;
s->nrx=l;
}
/* product */
for(j=0;j<r->npr;j++){
k=find_species(r->pr[j]);
if(k==-1){
printf("WARNING: species %s not found \n",r->re[j]);
continue;
}
s=X_spec+k;
l=s->nrx;
s->r[l]=i;
s->s[l]=r->spr[j]; /* change sign for reactant */
l++;
s->nrx=l;
}
}
}
/* the following code is probably suboptimal
it takes care of long names
*/
check_name_len(char *s)
{
char temp[9],x[5];
if(strlen(s)>9){
long_names[lnum].src=(char *)malloc(strlen(s)+1);
strcpy(long_names[lnum].src,s);
strncpy(x,s,4);
x[4]=0;
sprintf(long_names[lnum].rep,"%s.%d",x,lnum);
printf("long name: %s -> %s \n",long_names[lnum].src,long_names[lnum].rep);
lnum++;
}
}
/* replaces copies snew = sold with sfnd replaced by srep */
strrep(char *sold,char *snew,char *sfnd,char *srep)
{
int i=0,j=0,k,nf=strlen(sfnd),nr=strlen(srep);
int m=strlen(sold);
int l=0,lold=0;
while(1){
l=strfnd(sfnd,sold,lold);
if(l==-1){
for(k=lold;k<m;k++)
snew[i+k-lold]=sold[k];
snew[i+m-lold]=0;
return;
}
if(l>lold){
for(k=lold;k<l;k++){
snew[i]=sold[k];
i++;
}
}
/* okay it is there */
for(k=0;k<nr;k++){
snew[i]=srep[k];
i++;
}
lold=l+nf;
}
}
/* finds s1 in s2 starting at j0 */
int strfnd(char *s1,char *s2,int j0)
{
int r=-1,false;
int i=0,k=0,l;
int n1=strlen(s1),n2=strlen(s2);
if (n2<(n1+j0))return -1; /* cant happen */
while(1){
if(s2[j0+i]==s1[0]){
if((j0+i+n1)>n2)return -1; /* cant be included */
false=0;
l=i+j0;
for(k=0;k<n1;k++){
/* printf("i=%d,k=%d,%c %c \n",i,k,s2[j0+i],s1[k]); */
if(s2[j0+i]!=s1[k]){
false=1;
break;
}
i++;
}
if(false==0)return(l);
}
i++;
if(i>=n2)return -1;
}
}
fix_long_names(char *big,char *bigp)
{
int i=0;
char z[2048],zp[2048];
strcpy(z,big);
for(i=0;i<lnum;i++){
strrep(z,zp,long_names[i].src,long_names[i].rep);
strcpy(z,zp);
}
strcpy(bigp,z);
}
static int z_sort(sy1,sy2)
LONG_NAMES *sy1,*sy2;
{
if(strlen(sy1->src)>strlen(sy2->src))return -1;
return 1;
}
sort_long_names()
{
int i;
if(lnum<2)return; /* nothing to sort ! */
qsort(long_names,lnum,sizeof(LONG_NAMES),z_sort);
for(i=0;i<lnum;i++)
printf("%d: %s -> %s \n",i,long_names[i].src,long_names[i].rep);
}
write_ode_file(char *base)
{
char fname[256],tmp[2048];
char big[2048],bigp[2048];
FILE *fp;
RULE *r;
RXN *rx;
FUN_DEF *fn;
SPECIES *x;
EVENT *ev;
int i,j,k,na;
sprintf(fname,"%s.ode",base);
fp=fopen(fname,"w");
fprintf(fp,"# %s\n",fname);
fprintf(fp,"# Translated from %s by s2c \n",base);
fprintf(fp,"pow(x,y)=x^y\n");
fprintf(fp,"root(x,y)=y^(1/x)\n");
if(Nevent>0){
fprintf(fp,"gt(x,y)=x-y\n");
fprintf(fp,"lt(x,y)=y-x\n");
}
mark_rule_pars(); /* mark all parameters and species which
are actually rules so we dont double up */
sort_long_names(); /* reorder all names */
for(i=0;i<Nfuns;i++){
fn=funs+i;
na=fn->nargs;
sprintf(big,"%s(",fn->name);
for(j=0;j<na-1;j++){
sprintf(tmp,"%s,",fn->arg[j]);
strcat(big,tmp);
}
sprintf(tmp,"%s)=%s",fn->arg[na-1],fn->formula);
strcat(big,tmp);
fix_long_names(big,bigp);
fprintf(fp,"%s\n",bigp);
}
for(i=0;i<Nrule;i++){
r=rule+i;
if(!is_blank(r->v)){
sprintf(big,"%s=%s",r->v,r->f);
fix_long_names(big,bigp);
fprintf(fp,"%s\n",bigp);
}
/* dont print blank named rules, probably they are something else
like globals which we will find later
*/
}
for(i=0;i<Npar;i++){
if((par[i].fixed==-2)||(par[i].unique==-1))continue;
if(!is_blank(par[i].id))
sprintf(big,"par %s=%g",par[i].id,par[i].z);
else
sprintf(big,"par %s=%g",par[i].name,par[i].z);
fix_long_names(big,bigp);
fprintf(fp,"%s\n",bigp);
}
for(i=0;i<N_spec;i++){
x=X_spec+i;
if(x->rule==1)continue;
if(x->bc==1){
if(!is_blank(x->id))
sprintf(big,"%s=%g",x->id,x->x0);
else
sprintf(big,"%s=%g",x->name,x->x0);
}
else {
if(!is_blank(x->id))
sprintf(big,"init %s=%g",x->id,x->x0);
else
sprintf(big,"init %s=%g",x->name,x->x0);
}
fix_long_names(big,bigp);
fprintf(fp,"%s\n",bigp);
}
for(i=0;i<Nrxn;i++){
rx=rxn+i;
sprintf(big,"Rxn%d=%s",i+1,rx->formula);
fix_long_names(big,bigp);
fprintf(fp,"%s\n",bigp);
}
species_participation();
/* Finally write the damned equations ! */
for(i=0;i<N_spec;i++){
x=X_spec+i;
if(x->bc==1||x->rule==1)continue; /* dont do boundary conditions
or rules */
if(!is_blank(x->id))
sprintf(big,"d%s/dt=",x->id);
else
sprintf(big,"d%s/dt=",x->name);
for(j=0;j<x->nrx;j++){
k=x->r[j];
if(j>0){sprintf(tmp," + ");strcat(big,tmp);}
sprintf(tmp,"(%g)*Rxn%d",x->s[j],k+1);
strcat(big,tmp);
}
if(x->nrx==0){
sprintf(tmp,"0");
strcat(big,tmp);
}/* just to be Ok */
fix_long_names(big,bigp);
fprintf(fp,"%s\n",bigp);
}
/* last but not least, the globals */
for(i=0;i<Nevent;i++){
ev=event+i;
sprintf(big,"global 1 %s {",ev->ev);
na=ev->na;
for(j=0;j<na-1;j++){
sprintf(tmp,"%s;",ev->a[j]);
strcat(big,tmp);
}
sprintf(tmp,"%s}",ev->a[na-1]);
strcat(big,tmp);
fix_long_names(big,bigp);
fprintf(fp,"%s\n",bigp);
}
fprintf(fp,"done\n");
fclose(fp);
}
syntax highlighted by Code2HTML, v. 0.9.1