1 /* Steven Andrews, started 10/22/2001.
2  This is a library of functions for the Smoldyn program.
3  See documentation called SmoldynManual.pdf and SmoldynCodeDoc.pdf, and the Smoldyn
4  website, which is at www.smoldyn.org.
5  Copyright 2003-2016 by Steven Andrews.  This work is distributed under the terms
6  of the Gnu Lesser General Public License (LGPL). */
7 
8 #include <stdio.h>
9 #include <stdlib.h>
10 #include <math.h>
11 #include <string.h>
12 #include "string2.h"
13 #include "smoldyn.h"
14 #include "smoldynfuncs.h"
15 
16 #ifdef WINDOWS_BUILD
17 	#define DEVNULL "> NUL"
18 #else
19 	#define DEVNULL "> /dev/null"
20 #endif
21 
22 
23 /******************************************************************************/
24 /*********************************** BioNetGen ********************************/
25 /******************************************************************************/
26 
27 
28 
29 /******************************************************************************/
30 /****************************** Local declarations ****************************/
31 /******************************************************************************/
32 
33 // memory management
34 
35 bngptr bngalloc(bngptr bng,int maxparams,int maxmonomer,int maxbspecies,int maxbrxns);
36 void bngfree(bngptr bng);
37 bngssptr bngssalloc(bngssptr bngss,int maxbng);
38 void bngssfree(bngssptr bngss);
39 
40 // structure setup
41 
42 void bngsetcondition(bngssptr bngss,enum StructCond cond,int upgrade);
43 int bngenablebng(simptr sim,int maxbng);
44 bngptr bngaddbng(simptr sim,const char *bngname);
45 int bngparseparameter(bngptr bng,int index);
46 int bngaddparameter(bngptr bng,const char *name,const char *string);
47 int bngaddmonomer(bngptr bng,const char *name,enum MolecState ms);
48 
49 int bngmakeshortname(bngptr bng,int index,int totalmn,int hasmods);
50 enum MolecState bngmakedefaultstate(bngptr bng,int index,int totalmn);
51 double bngmakedifc(bngptr bng,int index,int totalmn);
52 
53 int bngaddspecies(bngptr bng,int bindex,const char *longname,const char *countstr);
54 int bngaddreaction(bngptr bng,int bindex,const char *reactants,const char *products,const char *rate);
55 int bngaddgroup(bngptr bng,int gindex,const char *gname,const char *specieslist);
56 bngptr bngreadstring(simptr sim,ParseFilePtr pfp,bngptr bng,const char *word,char *line2);
57 int bngupdateparams(simptr sim);
58 int bngupdatelists(simptr sim);
59 int bngupdate(simptr sim);
60 
61 
62 
63 /******************************************************************************/
64 /******************************* memory management ****************************/
65 /******************************************************************************/
66 
67 
68 /* bngallocsurfacedata */
bngallocsurfacedata(bngptr bng,int maxsurface)69 void bngallocsurfacedata(bngptr bng,int maxsurface) {
70 	int oldmaxsurface,s,i;
71 	enum SrfAction **newmonomeraction;
72 	surfactionptr **newmonomeractdetails;
73 	enum PanelFace face;
74 
75 	oldmaxsurface=bng->bngmaxsurface;
76 	for(i=0;i<bng->maxmonomer;i++) {
77 		if(!bng->monomeraction[i]) {
78 			CHECKMEM(bng->monomeraction[i]=(enum SrfAction**) calloc(maxsurface,sizeof(enum SrfAction*)));
79 			CHECKMEM(bng->monomeractdetails[i]=(surfactionptr **) calloc(maxsurface,sizeof(surfactionptr*)));
80 			for(s=0;s<maxsurface;s++) {
81 				bng->monomeraction[i][s]=NULL;
82 				bng->monomeractdetails[i][s]=NULL; }
83 			for(s=0;s<maxsurface;s++) {
84 				CHECKMEM(bng->monomeraction[i][s]=(enum SrfAction*) calloc(3,sizeof(enum SrfAction)));
85 				CHECKMEM(bng->monomeractdetails[i][s]=(surfactionptr*) calloc(3,sizeof(surfactionptr)));
86 				for(face=PFfront;face<=PFnone;face=(enum PanelFace)(face+1)) {
87 					bng->monomeraction[i][s][face]=SAtrans;
88 					bng->monomeractdetails[i][s][face]=NULL; }}}
89 
90 		else if(oldmaxsurface<maxsurface) {
91 			CHECKMEM(newmonomeraction=(enum SrfAction**) calloc(maxsurface,sizeof(enum SrfAction*)));
92 			CHECKMEM(newmonomeractdetails=(surfactionptr **) calloc(maxsurface,sizeof(surfactionptr*)));
93 			for(s=0;s<oldmaxsurface;s++) {
94 				newmonomeraction=bng->monomeraction[i];
95 				newmonomeractdetails=bng->monomeractdetails[i]; }
96 			for(s=oldmaxsurface;s<maxsurface;s++) {
97 				newmonomeraction[s]=NULL;
98 				newmonomeractdetails[s]=NULL; }
99 			for(s=oldmaxsurface;s<maxsurface;s++) {
100 				CHECKMEM(newmonomeraction[s]=(enum SrfAction*) calloc(3,sizeof(enum SrfAction)));
101 				CHECKMEM(newmonomeractdetails[s]=(surfactionptr*) calloc(3,sizeof(surfactionptr)));
102 				for(face=PFfront;face<=PFnone;face=(enum PanelFace)(face+1)) {
103 					newmonomeraction[s][face]=SAtrans;
104 					newmonomeractdetails[s][face]=NULL; }}
105 			free(bng->monomeraction[i]);
106 			free(bng->monomeractdetails[i]);
107 			bng->monomeraction[i]=newmonomeraction;
108 			bng->monomeractdetails[i]=newmonomeractdetails; }}
109 
110 		bng->bngmaxsurface=maxsurface;
111 	return;
112  failure:
113 	return; }
114 
115 
116 /* bngalloc */
bngalloc(bngptr bng,int maxparams,int maxmonomer,int maxbspecies,int maxbrxns)117 bngptr bngalloc(bngptr bng,int maxparams,int maxmonomer,int maxbspecies,int maxbrxns) {
118 	int i,newmax,oldmax;
119 	char **strlist;
120 	double *dbllist,**dblptrlist;
121   int *intlist,**intptrlist;
122   rxnptr *rxnlist;
123 	enum MolecState *mslist;
124 	enum SrfAction ***salist;
125 	surfactionptr ***sadlist;
126 
127 	newmax=oldmax=0;
128 	strlist=NULL;
129 	dbllist=NULL;
130   intlist=NULL;
131 
132 	if(!bng) {
133 		bng=(bngptr)malloc(sizeof(struct bngstruct));
134 		CHECKMEM(bng);
135     bng->bngss=NULL;
136 		bng->bngname=NULL;
137     bng->bngindex=0;
138 		bng->unirate=1;
139 		bng->birate=1;
140 		bng->maxparams=0;
141 		bng->nparams=0;
142 		bng->paramnames=NULL;
143     bng->paramstrings=NULL;
144 		bng->paramvalues=NULL;
145     bng->maxmonomer=0;
146     bng->nmonomer=0;
147     bng->monomernames=NULL;
148     bng->monomercount=NULL;
149 		bng->monomerdifc=NULL;
150 		bng->monomerdisplaysize=NULL;
151 		bng->monomercolor=NULL;
152 		bng->monomerstate=NULL;
153 		bng->bngmaxsurface=0;
154 		bng->monomeraction=NULL;
155 		bng->monomeractdetails=NULL;
156 		bng->maxbspecies=0;
157 		bng->nbspecies=0;
158 		bng->bsplongnames=NULL;
159 		bng->bspshortnames=NULL;
160 		bng->bspstate=NULL;
161 		bng->bspcountstr=NULL;
162 		bng->bspcount=NULL;
163 		bng->spindex=NULL;
164 		bng->maxbrxns=0;
165 		bng->nbrxns=0;
166 		bng->brxnreactstr=NULL;
167 		bng->brxnprodstr=NULL;
168     bng->brxnratestr=NULL;
169     bng->brxnreact=NULL;
170     bng->brxnprod=NULL;
171     bng->brxnorder=NULL;
172     bng->brxnnprod=NULL;
173     bng->brxn=NULL; }
174 
175 	if(maxparams>bng->maxparams) {					// enlarge params
176 		newmax=maxparams;
177 		oldmax=bng->maxparams;
178 		CHECKMEM(strlist=(char**)calloc(newmax,sizeof(char *)));
179 		for(i=0;i<oldmax;i++) strlist[i]=bng->paramnames[i];
180 		for(;i<newmax;i++) {CHECKMEM(strlist[i]=EmptyString());}
181 		free(bng->paramnames);
182 		bng->paramnames=strlist;
183 
184 		CHECKMEM(strlist=(char**)calloc(newmax,sizeof(char *)));
185 		for(i=0;i<oldmax;i++) strlist[i]=bng->paramstrings[i];
186 		for(;i<newmax;i++) {CHECKMEM(strlist[i]=EmptyString());}
187 		free(bng->paramstrings);
188 		bng->paramstrings=strlist;
189 
190 		CHECKMEM(dbllist=(double*)calloc(newmax,sizeof(double)));
191 		for(i=0;i<oldmax;i++) dbllist[i]=bng->paramvalues[i];
192 		for(;i<newmax;i++) dbllist[i]=0;
193 		free(bng->paramvalues);
194 		bng->paramvalues=dbllist;
195 
196     bng->maxparams=newmax; }
197 
198   if(maxmonomer>bng->maxmonomer) {          // enlarge monomers
199 		newmax=maxmonomer;
200 		oldmax=bng->maxmonomer;
201 		CHECKMEM(strlist=(char**)calloc(newmax,sizeof(char *)));
202 		for(i=0;i<oldmax;i++) strlist[i]=bng->monomernames[i];
203 		for(;i<newmax;i++) {CHECKMEM(strlist[i]=EmptyString());}
204 		free(bng->monomernames);
205 		bng->monomernames=strlist;
206 
207     CHECKMEM(intlist=(int*)calloc(newmax,sizeof(int)));
208 		for(i=0;i<oldmax;i++) intlist[i]=bng->monomercount[i];
209 		for(;i<newmax;i++) intlist[i]=0;
210 		free(bng->monomercount);
211 		bng->monomercount=intlist;
212 
213 		CHECKMEM(dbllist=(double*)calloc(newmax,sizeof(double)));
214 		for(i=0;i<oldmax;i++) dbllist[i]=bng->monomerdifc[i];
215 		for(;i<newmax;i++) dbllist[i]=0;
216 		free(bng->monomerdifc);
217 		bng->monomerdifc=dbllist;
218 
219 		CHECKMEM(dbllist=(double*)calloc(newmax,sizeof(double)));
220 		for(i=0;i<oldmax;i++) dbllist[i]=bng->monomerdisplaysize[i];
221 		for(;i<newmax;i++) dbllist[i]=0;
222 		free(bng->monomerdisplaysize);
223 		bng->monomerdisplaysize=dbllist;
224 
225 		CHECKMEM(dblptrlist=(double**)calloc(newmax,sizeof(double *)));
226 		for(i=0;i<oldmax;i++) dblptrlist[i]=bng->monomercolor[i];
227 		for(;i<newmax;i++) {
228 			CHECKMEM(dblptrlist[i]=(double*)calloc(3,sizeof(double)));
229 			dblptrlist[i][0]=0;
230 			dblptrlist[i][1]=0;
231 			dblptrlist[i][2]=0; }
232 		free(bng->monomercolor);
233 		bng->monomercolor=dblptrlist;
234 
235 		CHECKMEM(mslist=(enum MolecState *)calloc(newmax,sizeof(enum MolecState)));
236 		for(i=0;i<oldmax;i++) mslist[i]=bng->monomerstate[i];
237 		for(;i<newmax;i++) mslist[i]=MSsoln;
238 		free(bng->monomerstate);
239 		bng->monomerstate=mslist;
240 
241 		CHECKMEM(salist=(enum SrfAction ***)calloc(newmax,sizeof(enum SrfAction **)));
242 		for(i=0;i<oldmax;i++) salist[i]=bng->monomeraction[i];
243 		for(;i<newmax;i++) salist[i]=NULL;
244 		free(bng->monomeraction);
245 		bng->monomeraction=salist;
246 
247 		CHECKMEM(sadlist=(surfactionptr ***)calloc(newmax,sizeof(surfactionptr **)));
248 		for(i=0;i<oldmax;i++) sadlist[i]=bng->monomeractdetails[i];
249 		for(;i<newmax;i++) sadlist[i]=NULL;
250 		free(bng->monomeractdetails);
251 		bng->monomeractdetails=sadlist;
252 
253 		bng->maxmonomer=newmax; }
254 
255 	if(maxbspecies>bng->maxbspecies) {				// enlarge bspecies
256 		newmax=maxbspecies;
257 		oldmax=bng->maxbspecies;
258 		CHECKMEM(strlist=(char**)calloc(newmax,sizeof(char *)));
259 		for(i=0;i<oldmax;i++) strlist[i]=bng->bsplongnames[i];
260 		for(;i<newmax;i++) {CHECKMEM(strlist[i]=EmptyString());}
261 		free(bng->bsplongnames);
262 		bng->bsplongnames=strlist;
263 
264     CHECKMEM(strlist=(char**)calloc(newmax,sizeof(char *)));
265 		for(i=0;i<oldmax;i++) strlist[i]=bng->bspshortnames[i];
266 		for(;i<newmax;i++) {CHECKMEM(strlist[i]=EmptyString());}
267 		free(bng->bspshortnames);
268 		bng->bspshortnames=strlist;
269 
270 		CHECKMEM(mslist=(enum MolecState *)calloc(newmax,sizeof(enum MolecState)));
271 		for(i=0;i<oldmax;i++) mslist[i]=bng->bspstate[i];
272 		for(;i<newmax;i++) mslist[i]=MSsoln;
273 		free(bng->bspstate);
274 		bng->bspstate=mslist;
275 
276     CHECKMEM(strlist=(char**)calloc(newmax,sizeof(char *)));
277 		for(i=0;i<oldmax;i++) strlist[i]=bng->bspcountstr[i];
278 		for(;i<newmax;i++) {CHECKMEM(strlist[i]=EmptyString());}
279 		free(bng->bspcountstr);
280 		bng->bspcountstr=strlist;
281 
282 		CHECKMEM(dbllist=(double*)calloc(newmax,sizeof(double)));
283 		for(i=0;i<oldmax;i++) dbllist[i]=bng->bspcount[i];
284 		for(;i<newmax;i++) dbllist[i]=0;
285 		free(bng->bspcount);
286 		bng->bspcount=dbllist;
287 
288     CHECKMEM(intlist=(int*)calloc(newmax,sizeof(int)));
289 		for(i=0;i<oldmax;i++) intlist[i]=bng->spindex[i];
290 		for(;i<newmax;i++) intlist[i]=0;
291 		free(bng->spindex);
292 		bng->spindex=intlist;
293 
294     bng->maxbspecies=newmax; }
295 
296 	if(maxbrxns>bng->maxbrxns) {				// enlarge brxns
297 		newmax=maxbrxns;
298 		oldmax=bng->maxbrxns;
299 		CHECKMEM(strlist=(char**)calloc(newmax,sizeof(char *)));
300 		for(i=0;i<oldmax;i++) strlist[i]=bng->brxnreactstr[i];
301 		for(;i<newmax;i++) {CHECKMEM(strlist[i]=EmptyString());}
302 		free(bng->brxnreactstr);
303 		bng->brxnreactstr=strlist;
304 
305 		CHECKMEM(strlist=(char**)calloc(newmax,sizeof(char *)));
306 		for(i=0;i<oldmax;i++) strlist[i]=bng->brxnprodstr[i];
307 		for(;i<newmax;i++) {CHECKMEM(strlist[i]=EmptyString());}
308 		free(bng->brxnprodstr);
309 		bng->brxnprodstr=strlist;
310 
311 		CHECKMEM(strlist=(char**)calloc(newmax,sizeof(char *)));
312 		for(i=0;i<oldmax;i++) strlist[i]=bng->brxnratestr[i];
313 		for(;i<newmax;i++) {CHECKMEM(strlist[i]=EmptyString());}
314 		free(bng->brxnratestr);
315 		bng->brxnratestr=strlist;
316 
317     CHECKMEM(intptrlist=(int**)calloc(newmax,sizeof(int *)));
318 		for(i=0;i<oldmax;i++) intptrlist[i]=bng->brxnreact[i];
319 		for(;i<newmax;i++) {
320       CHECKMEM(intptrlist[i]=(int*)calloc(2,sizeof(int)));
321       intptrlist[i][0]=0;
322       intptrlist[i][1]=0; }
323 		free(bng->brxnreact);
324 		bng->brxnreact=intptrlist;
325 
326     CHECKMEM(intptrlist=(int**)calloc(newmax,sizeof(int *)));
327 		for(i=0;i<oldmax;i++) intptrlist[i]=bng->brxnprod[i];
328 		for(;i<newmax;i++) {
329       CHECKMEM(intptrlist[i]=(int*)calloc(2,sizeof(int)));
330       intptrlist[i][0]=0;
331       intptrlist[i][1]=0; }
332 		free(bng->brxnprod);
333 		bng->brxnprod=intptrlist;
334 
335     CHECKMEM(intlist=(int*)calloc(newmax,sizeof(int)));
336 		for(i=0;i<oldmax;i++) intlist[i]=bng->brxnorder[i];
337 		for(;i<newmax;i++) intlist[i]=0;
338 		free(bng->brxnorder);
339 		bng->brxnorder=intlist;
340 
341     CHECKMEM(intlist=(int*)calloc(newmax,sizeof(int)));
342 		for(i=0;i<oldmax;i++) intlist[i]=bng->brxnnprod[i];
343 		for(;i<newmax;i++) intlist[i]=0;
344 		free(bng->brxnnprod);
345 		bng->brxnnprod=intlist;
346 
347     CHECKMEM(rxnlist=(rxnptr*)calloc(newmax,sizeof(rxnptr)));
348 		for(i=0;i<oldmax;i++) rxnlist[i]=bng->brxn[i];
349 		for(;i<newmax;i++) rxnlist[i]=NULL;
350 		free(bng->brxn);
351 		bng->brxn=rxnlist;
352 
353     bng->maxbrxns=newmax; }
354 
355 		return bng;
356 failure:
357 	bngfree(bng);
358 	simLog(NULL,10,"Unable to allocate memory in bngalloc");
359 	return NULL; }
360 
361 
362 /* bngfree */
bngfree(bngptr bng)363 void bngfree(bngptr bng) {
364 	int i,s;
365 
366 	if(!bng) return;
367 
368 	for(i=0;i<bng->maxbrxns;i++) {
369 		free(bng->brxnreactstr[i]);
370 		free(bng->brxnprodstr[i]);
371     free(bng->brxnratestr[i]);
372     free(bng->brxnreact[i]);
373     free(bng->brxnprod[i]); }
374 	free(bng->brxnreactstr);
375 	free(bng->brxnprodstr);
376   free(bng->brxnreact);
377   free(bng->brxnprod);
378   free(bng->brxnorder);
379   free(bng->brxnnprod);
380   free(bng->brxn);
381 
382 	for(i=0;i<bng->maxbspecies;i++) {
383 		free(bng->bsplongnames[i]);
384 		free(bng->bspshortnames[i]);
385 		free(bng->bspcountstr[i]); }
386 	free(bng->bsplongnames);
387 	free(bng->bspshortnames);
388 	free(bng->bspstate);
389 	free(bng->bspcountstr);
390 	free(bng->bspcount);
391 	free(bng->spindex);
392 
393   for(i=0;i<bng->maxmonomer;i++) {
394     free(bng->monomernames[i]); }
395   free(bng->monomernames);
396   free(bng->monomercount);
397 	free(bng->monomerdifc);
398 	free(bng->monomerdisplaysize);
399 	for(i=0;i<bng->maxmonomer;i++)
400 		free(bng->monomercolor[i]);
401 	free(bng->monomercolor);
402 	free(bng->monomerstate);
403 	for(i=0;i<bng->maxmonomer;i++) {
404 		if(bng->monomeraction[i])
405 			for(s=0;s<bng->bngmaxsurface;s++) {
406 				free(bng->monomeraction[i][s]);
407 				free(bng->monomeractdetails[i][s]); }
408 		free(bng->monomeraction[i]);
409 		free(bng->monomeractdetails[i]); }
410 	free(bng->monomeraction);
411 	free(bng->monomeractdetails);
412 
413 	for(i=0;i<bng->maxparams;i++) {
414 		free(bng->paramnames[i]);
415 		free(bng->paramstrings[i]); }
416 	free(bng->paramnames);
417 	free(bng->paramstrings);
418 	free(bng->paramvalues);
419 
420 	free(bng);
421 	return; }
422 
423 
424 /* bngssalloc */
bngssalloc(bngssptr bngss,int maxbng)425 bngssptr bngssalloc(bngssptr bngss,int maxbng) {
426 	int i,newmax,oldmax;
427 	char **strlist;
428 	bngptr *newbnglist;
429 
430 	if(!bngss) {
431 		CHECKMEM(bngss=(bngssptr) malloc(sizeof(struct bngsuperstruct)));
432 		bngss->condition=SCinit;
433 		bngss->sim=NULL;
434 		bngss->BNG2path=NULL;
435 		bngss->maxbng=0;
436 		bngss->nbng=0;
437 		bngss->bngnames=NULL;
438 		bngss->bnglist=NULL; }
439 
440 	if(!bngss->BNG2path) {
441 		CHECKMEM(bngss->BNG2path=EmptyString());
442 		strcpy(bngss->BNG2path,BNG2_PATH); }
443 
444 	if(maxbng>bngss->maxbng) {
445 		newmax=maxbng;
446 		oldmax=bngss->maxbng;
447     CHECKMEM(strlist=(char**)calloc(newmax,sizeof(char *)));
448 		for(i=0;i<oldmax;i++) strlist[i]=bngss->bngnames[i];
449 		for(;i<newmax;i++) {CHECKMEM(strlist[i]=EmptyString());}
450 		free(bngss->bngnames);
451 		bngss->bngnames=strlist;
452 
453     CHECKMEM(newbnglist=(bngptr*)calloc(newmax,sizeof(bngptr)));
454 		for(i=0;i<oldmax;i++) newbnglist[i]=bngss->bnglist[i];
455 		for(;i<newmax;i++) {
456 			CHECKMEM(newbnglist[i]=bngalloc(NULL,1,1,1,1));
457 			newbnglist[i]->bngss=bngss;
458 			newbnglist[i]->bngname=bngss->bngnames[i];
459       newbnglist[i]->bngindex=i; }
460 		free(bngss->bnglist);
461 		bngss->bnglist=newbnglist;
462 
463     bngss->maxbng=newmax; }
464 	return bngss;
465 
466 failure:
467 	bngssfree(bngss);
468 	simLog(NULL,10,"Unable to allocate memory in bngssalloc");
469 	return NULL; }
470 
471 
472 /* bngssfree */
bngssfree(bngssptr bngss)473 void bngssfree(bngssptr bngss) {
474 	int i;
475 
476 	if(!bngss) return;
477 	for(i=0;i<bngss->maxbng;i++) {
478 		bngfree(bngss->bnglist[i]);
479 		free(bngss->bngnames[i]); }
480 	free(bngss->bnglist);
481 	free(bngss->bngnames);
482 	free(bngss->BNG2path);
483 	free(bngss);
484 	return; }
485 
486 
487 /******************************************************************************/
488 /***************************** data structure output **************************/
489 /******************************************************************************/
490 
491 
492 /* bngoutput */
bngoutput(simptr sim)493 void bngoutput(simptr sim) {
494 	bngssptr bngss;
495 	bngptr bng;
496 	int b,i,s;
497 	char string[STRCHAR];
498 
499 	bngss=sim->bngss;
500 	if(!bngss) return;
501 	simLog(sim,2,"BioNetGen parameters\n");
502 	simLog(sim,2," BNG2.pl path: %s\n",bngss->BNG2path);
503 	simLog(sim,1," BNG allocated: %i,",bngss->maxbng);
504 	simLog(sim,2," BNG defined: %i\n",bngss->nbng);
505 	for(b=0;b<bngss->nbng;b++) {
506 		bng=bngss->bnglist[b];
507 		simLog(sim,2," BNG: %s\n",bng->bngname);
508 		if(bng->unirate!=1 || bng->birate!=1)
509 			simLog(sim,2,"  rate multipliers: unimolecular: %g, bimolecular: %g\n",bng->unirate,bng->birate);
510 		simLog(sim,1,"  parameters allocated: %i,",bng->maxparams);
511 		simLog(sim,2,"  parameters defined: %i\n",bng->nparams);
512 		for(i=0;i<bng->nparams;i++) {
513 			simLog(sim,2,"   %i %s %g\n",i,bng->paramnames[i],bng->paramvalues[i]); }
514 		simLog(sim,1,"  monomers allocated: %i,",bng->maxmonomer);
515 		simLog(sim,2,"  monomers defined: %i\n",bng->nmonomer);
516 		for(i=0;i<bng->nmonomer;i++) {
517 			simLog(sim,2,"   %s: default state: %s, diffusion coeff.: %g\n",bng->monomernames[i],molms2string(bng->monomerstate[i],string),bng->monomerdifc[i]);
518 			simLog(sim,2,"    display size: %g, color: %g %g %g\n",bng->monomerdisplaysize[i],bng->monomercolor[i][0],bng->monomercolor[i][1],bng->monomercolor[i][2]);
519 			for(s=0;s<bng->bngmaxsurface;s++) {
520 				simLog(sim,2,"    for surface %s: %s at front",sim->srfss->srflist[s]->sname,surfact2string(bng->monomeraction[i][s][PFfront],string));
521 				simLog(sim,2,", %s at back\n",surfact2string(bng->monomeraction[i][s][PFback],string)); }}
522 		simLog(sim,1,"  species allocated: %i,",bng->maxbspecies);
523 		simLog(sim,2,"  species defined: %i\n",bng->nbspecies-1);
524 		for(i=0;i<bng->nbspecies;i++) {
525 			if(bng->spindex[i]>0)
526 				simLog(sim,2,"   %i %s (%s), count: %g, longname: %s\n",i,bng->bspshortnames[i],molms2string(bng->bspstate[i],string),bng->bspcount[i],bng->bsplongnames[i]); }
527 		simLog(sim,1,"  reactions allocated: %i,",bng->maxbrxns);
528 		simLog(sim,2,"  reactions defined: %i\n",bng->nbrxns-1);
529 		for(i=0;i<bng->nbrxns;i++)
530 			if(bng->brxn[i]) {
531 				simLog(sim,2,"   %i",i);
532 				if(bng->brxnorder[i]>=1) simLog(sim,2," %s",bng->bspshortnames[bng->brxnreact[i][0]]);
533 				if(bng->brxnorder[i]==2) simLog(sim,2," + %s",bng->bspshortnames[bng->brxnreact[i][1]]);
534 				simLog(sim,2," ->");
535 				if(bng->brxnnprod[i]>=1) simLog(sim,2," %s",bng->bspshortnames[bng->brxnprod[i][0]]);
536 				if(bng->brxnnprod[i]==2) simLog(sim,2," + %s",bng->bspshortnames[bng->brxnprod[i][1]]);
537 				simLog(sim,2,"  rate: %g",bng->brxn[i]->rate);
538 				simLog(sim,2,"\n"); }}
539 	simLog(sim,2,"\n");
540 	return; }
541 
542 
543 /* checkbngparams */
checkbngparams(simptr sim,int * warnptr)544 int checkbngparams(simptr sim,int *warnptr) {
545 	int error,warn,b,i;
546 	bngssptr bngss;
547 	bngptr bng;
548 	char string[STRCHAR];
549 
550 	error=warn=0;
551 	bngss=sim->bngss;
552 	if(!bngss) {
553 		if(warnptr) *warnptr=warn;
554 		return 0; }
555 
556 	if(bngss->condition!=SCok) {
557 		warn++;
558 		simLog(sim,7," WARNING: bng structure condition is %s\n",simsc2string(bngss->condition,string)); }
559 
560 	for(b=0;b<bngss->nbng;b++) {
561 		bng=bngss->bnglist[b];
562 		for(i=0;i<bng->nbspecies;i++) {
563 			if(bng->bspcount[i]>0 && bng->bspcount[i]<1) simLog(sim,7," WARNING: count for %s is very low\n",bng->bspshortnames[i]); }}
564 
565 	if(warnptr) *warnptr=warn;
566 	return error; }
567 
568 
569 /******************************************************************************/
570 /**************************** structure set up - bng **************************/
571 /******************************************************************************/
572 
573 
574 /* bngsetcondition */
bngsetcondition(bngssptr bngss,enum StructCond cond,int upgrade)575 void bngsetcondition(bngssptr bngss,enum StructCond cond,int upgrade) {
576 	if(!bngss) return;
577 	if(upgrade==0 && bngss->condition>cond) bngss->condition=cond;
578 	else if(upgrade==1 && bngss->condition<cond) bngss->condition=cond;
579 	else if(upgrade==2) bngss->condition=cond;
580 	if(bngss->sim && bngss->condition<bngss->sim->condition) {
581 		cond=bngss->condition;
582 		simsetcondition(bngss->sim,cond==SCinit?SClists:cond,0); }
583 	return; }
584 
585 
586 /* bngenablebng */
bngenablebng(simptr sim,int maxbng)587 int bngenablebng(simptr sim,int maxbng) {
588 	bngssptr bngss;
589 
590 	if(sim->bngss)									// check for redundant function call
591 		if(maxbng==-1 || sim->bngss->maxbng==maxbng)
592 			return 0;
593 	bngss=bngssalloc(sim->bngss,maxbng<0?1:maxbng);
594 	if(!bngss) return 1;
595 	sim->bngss=bngss;
596 	bngss->sim=sim;
597 	bngsetcondition(sim->bngss,SClists,0);
598 	return 0; }
599 
600 
601 /* bngaddbng */
bngaddbng(simptr sim,const char * bngname)602 bngptr bngaddbng(simptr sim,const char *bngname) {
603 	int er,i;
604 	bngssptr bngss;
605 	bngptr bng;
606 
607 	if(!sim->bngss) {
608 	 er=bngenablebng(sim,-1);
609 	 if(er) return NULL; }
610 	bngss=sim->bngss;
611 
612 	i=stringfind(bngss->bngnames,bngss->nbng,bngname);
613 	if(i<0) {
614 		if(bngss->nbng==bngss->maxbng) {
615 			er=bngenablebng(sim,bngss->nbng*2+1);
616 			if(er) return NULL; }
617 		i=bngss->nbng++;
618 		strncpy(bngss->bngnames[i],bngname,STRCHAR-1);
619 		bngss->bngnames[i][STRCHAR-1]='\0';
620 		bng=bngss->bnglist[i]; }
621 	else
622 		bng=bngss->bnglist[i];
623 
624 	bngsetcondition(bngss,SClists,0);
625 	return bng; }
626 
627 
628 /* bngsetparam */
bngsetparam(bngptr bng,char * parameter,double amount)629 int bngsetparam(bngptr bng,char *parameter,double amount) {
630 	if(!strcmp(parameter,"unimolecular_rate")) {
631 		if(amount<0) return 2;
632 		bng->unirate=amount; }
633 	else if(!strcmp(parameter,"bimolecular_rate")) {
634 		if(amount<0) return 2;
635 		bng->birate=amount; }
636 	else
637 		return 1;
638 	return 0; }
639 
640 
641 /* bngsetBNG2path */
bngsetBNG2path(bngptr bng,char * path)642 int bngsetBNG2path(bngptr bng,char *path) {
643 	strcpy(bng->bngss->BNG2path,path);
644 	return 0; }
645 
646 
647 /******************************************************************************/
648 /************************ structure set up - parameters ***********************/
649 /******************************************************************************/
650 
651 
652 /* bngparseparameter */
bngparseparameter(bngptr bng,int index)653 int bngparseparameter(bngptr bng,int index) {
654   int er;
655   double value;
656 
657   er=0;
658   if(bng->paramstrings[index]) {
659     value=strmatheval(bng->paramstrings[index],bng->paramnames,bng->paramvalues,bng->nparams);
660     er=strmatherror(NULL,1);
661     bng->paramvalues[index]=value; }
662   return (er==0)?0:1; }
663 
664 
665 /* bngaddparameter */
bngaddparameter(bngptr bng,const char * name,const char * string)666 int bngaddparameter(bngptr bng,const char *name,const char *string) {
667 	int i,er;
668 
669 	if(bng->nparams==bng->maxparams) {		// enlarge list if needed
670 		bng=bngalloc(bng,bng->maxparams*2+1,0,0,0);
671 		if(!bng) return -1; }
672 
673   i=stringfind(bng->paramnames,bng->nparams,name);
674   if(i<0) {
675     i=bng->nparams;
676     bng->nparams++;
677     strcpy(bng->paramnames[i],name); }
678   if(string) strcpy(bng->paramstrings[i],string);
679   else bng->paramstrings[i][0]='\0';
680   er=bngparseparameter(bng,i);
681   return (er==0)?i:-2; }
682 
683 
684 /******************************************************************************/
685 /************************** structure set up - monomers ***********************/
686 /******************************************************************************/
687 
688 
689 /* bngaddmonomer */
bngaddmonomer(bngptr bng,const char * name,enum MolecState ms)690 int bngaddmonomer(bngptr bng,const char *name,enum MolecState ms) {
691 	simptr sim;
692   int i,i2,j,n,newcount,len,s;
693 	char *newname,*spname;
694 	double newdifc,newdisplaysize,*newcolor;
695 	enum MolecState newstate;
696 	surfaceptr srf;
697 	enum SrfAction **newaction;
698 	surfactionptr **newactdetails;
699 	enum PanelFace face;
700 
701 	sim=bng->bngss->sim;
702   i=stringfind(bng->monomernames,bng->nmonomer,name);	// see if it's already in there
703 
704   if(i<0) {
705 		if(!strokname(name))
706 			return -2;
707 		if(ms==MSbsoln) ms=MSsoln;
708 
709 		if(bng->nmonomer==bng->maxmonomer) {		// enlarge list if needed
710 			bng=bngalloc(bng,0,bng->maxmonomer*2+1,0,0);
711 			if(!bng) return -1; }
712 
713 		bng->nmonomer++;
714 		n=bng->nmonomer;
715 		strcpy(bng->monomernames[n-1],name);					// put the new monomer at the end
716 
717 		j=stringfind(sim->mols->spname,sim->mols->nspecies,name);	// look for monomer in Smoldyn species list
718 		if(j<=0) {																		// not there, so try Smoldyn species without dots
719 			for(j=1;j<sim->mols->nspecies;j++) {
720 				spname=sim->mols->spname[j];
721 				if(strsymbolcount(spname,'.')==2 && (len=strlen(name))==strchr(spname,'.')-spname && !strncmp(name,spname,len)) break; }
722 			if(j==sim->mols->nspecies) j=-1; }
723 		if(j>0) {																			// found Smoldyn species so add in data
724 			bng->monomerdifc[n-1]=sim->mols->difc[j][ms];
725 			bng->monomerdisplaysize[n-1]=sim->mols->display[j][ms];
726 			bng->monomercolor[n-1][0]=sim->mols->color[j][ms][0];
727 			bng->monomercolor[n-1][1]=sim->mols->color[j][ms][1];
728 			bng->monomercolor[n-1][2]=sim->mols->color[j][ms][2];
729 			if(sim->srfss) {														// assign surface interaction data
730 				if(bng->bngmaxsurface<sim->srfss->nsrf || !bng->monomeraction[n-1]) {
731 					bngallocsurfacedata(bng,sim->srfss->nsrf);
732 					if(bng->bngmaxsurface<sim->srfss->nsrf) return -1; }
733 				for(s=0;s<sim->srfss->nsrf;s++) {
734 					srf=sim->srfss->srflist[s];
735 					for(face=PFfront;face<=PFnone;face=(enum PanelFace)(face+1)) {
736 						bng->monomeraction[n-1][s][face]=srf->action[j][ms][face];
737 						bng->monomeractdetails[n-1][s][face]=srf->actdetails[j][ms][face]; }}}}
738 
739 		newname=bng->monomernames[n-1];								// create extra copy of the new stuff
740 		newcount=bng->monomercount[n-1];
741 		newdifc=bng->monomerdifc[n-1];
742 		newdisplaysize=bng->monomerdisplaysize[n-1];
743 		newcolor=bng->monomercolor[n-1];
744 		newstate=bng->monomerstate[n-1];
745 		newaction=bng->monomeraction[n-1];
746 		newactdetails=bng->monomeractdetails[n-1];
747 
748 		for(i=0;i<n-1;i++)														// sort new monomer to where it belongs
749 			if(strcmp(newname,bng->monomernames[i])<0) {
750 				for(i2=n-2;i2>=i;i2--) {
751 					bng->monomernames[i2+1]=bng->monomernames[i2];
752 					bng->monomercount[i2+1]=bng->monomercount[i2];
753 					bng->monomerdifc[i2+1]=bng->monomerdifc[i2];
754 					bng->monomerdisplaysize[i2+1]=bng->monomerdisplaysize[i2];
755 					bng->monomercolor[i2+1]=bng->monomercolor[i2];
756 					bng->monomerstate[i2+1]=bng->monomerstate[i2];
757 					bng->monomeraction[i2+1]=bng->monomeraction[i2];
758 					bng->monomeractdetails[i2+1]=bng->monomeractdetails[i2]; }
759 				bng->monomernames[i]=newname;
760 				bng->monomercount[i]=newcount;
761 				bng->monomerdifc[i]=newdifc;
762 				bng->monomerdisplaysize[i]=newdisplaysize;
763 				bng->monomercolor[i]=newcolor;
764 				bng->monomerstate[i]=newstate;
765 				bng->monomeraction[i]=newaction;
766 				bng->monomeractdetails[i]=newactdetails;
767 				break; }}
768 
769   return i; }
770 
771 
772 /* bngsetmonomerdifc */
bngsetmonomerdifc(bngptr bng,char * name,double difc)773 int bngsetmonomerdifc(bngptr bng,char *name,double difc) {
774 	int i;
775 
776 	if(!strcmp(name,"all")) {
777 		for(i=0;i<bng->nmonomer;i++)
778 			bng->monomerdifc[i]=difc; }
779 	else {
780 		i=bngaddmonomer(bng,name,MSsoln);
781 		if(i<0) return i;
782 		bng->monomerdifc[i]=difc; }
783 
784 	return 0; }
785 
786 
787 /* bngsetmonomerdisplaysize */
bngsetmonomerdisplaysize(bngptr bng,char * name,double displaysize)788 int bngsetmonomerdisplaysize(bngptr bng,char *name,double displaysize) {
789 	int i;
790 
791 	if(!strcmp(name,"all")) {
792 		for(i=0;i<bng->nmonomer;i++)
793 			bng->monomerdisplaysize[i]=displaysize; }
794 	else {
795 		i=bngaddmonomer(bng,name,MSsoln);
796 		if(i<0) return i;
797 		bng->monomerdisplaysize[i]=displaysize; }
798 
799 	return 0; }
800 
801 
802 /* bngsetmonomercolor */
bngsetmonomercolor(bngptr bng,char * name,double * color)803 int bngsetmonomercolor(bngptr bng,char *name,double *color) {
804 	int i;
805 
806 	if(!strcmp(name,"all")) {
807 		for(i=0;i<bng->nmonomer;i++) {
808 			bng->monomercolor[i][0]=color[0];
809 			bng->monomercolor[i][1]=color[1];
810 			bng->monomercolor[i][2]=color[2]; }}
811 	else {
812 		i=bngaddmonomer(bng,name,MSsoln);
813 		if(i<0) return i;
814 		bng->monomercolor[i][0]=color[0];
815 		bng->monomercolor[i][1]=color[1];
816 		bng->monomercolor[i][2]=color[2]; }
817 
818 	return 0; }
819 
820 
821 /* bngsetmonomerstate */
bngsetmonomerstate(bngptr bng,char * name,enum MolecState ms)822 int bngsetmonomerstate(bngptr bng,char *name,enum MolecState ms) {
823 	int i;
824 
825 	if(!strcmp(name,"all")) {
826 		for(i=0;i<bng->nmonomer;i++)
827 			bng->monomerstate[i]=ms; }
828 	else {
829 		i=bngaddmonomer(bng,name,ms);
830 		if(i<0) return i;
831 		bng->monomerstate[i]=ms; }
832 
833 	return 0; }
834 
835 
836 /******************************************************************************/
837 /************************** structure set up - species ************************/
838 /******************************************************************************/
839 
840 
841 /* bngmakeshortname */
bngmakeshortname(bngptr bng,int index,int totalmn,int hasmods)842 int bngmakeshortname(bngptr bng,int index,int totalmn,int hasmods) {
843 	char *shortname,string[STRCHAR],*lastdot,*cmpname;
844 	int length,mn,i1,i2,snlength,cmplength;
845 
846 	shortname=bng->bspshortnames[index];				// generate the short name root, diffusion coefficient, and default state
847 	shortname[0]='\0';
848 	length=STRCHAR-20;													// length remaining in shortname string
849 
850 	if(totalmn==1 && hasmods==0) {							// just 1 monomer and no modifications possible
851 		for(mn=0;mn<bng->nmonomer;mn++)
852 			if(bng->monomercount[mn]>0) {
853 				strcpy(shortname,bng->monomernames[mn]);
854 				mn=bng->nmonomer; }}
855 
856 	else {																			// multiple monomers or 1 monomer and possible modifications
857 		for(mn=0;mn<bng->nmonomer && length>0;mn++)
858 			if(bng->monomercount[mn]>0) {
859 				snprintf(string,STRCHAR,"%s.%i.",bng->monomernames[mn],bng->monomercount[mn]);
860 				string[length-1]='\0';
861 				strcat(shortname,string);
862 				length-=strlen(string); }
863 
864 		i2=0;																			// append the short name isomer number
865 		snlength=strlen(shortname);
866 		for(i1=0;i1<index;i1++) {
867 			cmpname=bng->bspshortnames[i1];
868 			lastdot=strrchr(cmpname,'.');
869 			if(lastdot) {
870 				cmplength=lastdot-cmpname;
871 				if(!strncmp(shortname,cmpname,cmplength>snlength?cmplength:snlength)) i2++; }}
872 		snprintf(string,STRCHAR,"%i",i2);
873 		strcat(shortname,string); }
874 
875 	return 0; }
876 
877 
878 /* bngmakedefaultstate */
bngmakedefaultstate(bngptr bng,int index,int totalmn)879 enum MolecState bngmakedefaultstate(bngptr bng,int index,int totalmn) {
880 	enum MolecState ms,trialms;
881 	int mn;
882 
883 	ms=MSsoln;
884 
885 	if(totalmn==1) {      // just 1 monomer
886 		for(mn=0;mn<bng->nmonomer;mn++)
887 			if(bng->monomercount[mn]>0) {
888 				ms=bng->monomerstate[mn];
889 				mn=bng->nmonomer; }}
890 
891 	else {                // multiple monomers
892 		for(mn=0;mn<bng->nmonomer;mn++)
893 			if(bng->monomercount[mn]>0) {
894 				trialms=bng->monomerstate[mn];
895 				if(ms==MSbsoln) {
896 					if(trialms!=MSsoln) ms=trialms; }
897 				else if(trialms==MSbsoln) {
898 					if(ms==MSsoln) ms=trialms; }
899 				else if(trialms>ms) ms=trialms; }}
900 
901 	bng->bspstate[index]=ms;
902 
903 	return ms; }
904 
905 
906 /* bngmakedifc */
bngmakedifc(bngptr bng,int index,int totalmn)907 double bngmakedifc(bngptr bng,int index,int totalmn) {
908 	double difc;
909 	int mn,j;
910 	simptr sim;
911 	enum MolecState ms;
912 
913 	difc=-1;
914 	sim=bng->bngss->sim;
915 	j=stringfind(sim->mols->spname,sim->mols->nspecies,bng->bspshortnames[index]);	// look for species in Smoldyn species list
916 
917 	if(j>0) {
918 		ms=bng->bspstate[index];
919 		if(ms==MSbsoln) ms=MSsoln;
920 		difc=sim->mols->difc[j][ms]; }
921 	else if(totalmn==1) {      // just 1 monomer
922 		for(mn=0;mn<bng->nmonomer;mn++)
923 			if(bng->monomercount[mn]>0) {
924 				difc=bng->monomerdifc[mn];
925 				mn=bng->nmonomer; }}
926 	else {                // multiple monomers
927 		for(mn=0;mn<bng->nmonomer;mn++)
928 			if(bng->monomercount[mn]>0) {
929 				if(difc==-1 && bng->monomerdifc[mn]==0) difc=0;
930 				else if(difc==-1) difc=bng->monomercount[mn]*pow(bng->monomerdifc[mn],-3.0);
931 				else if(bng->monomerdifc[mn]==0) difc=0;
932 				else difc+=bng->monomercount[mn]*pow(bng->monomerdifc[mn],-3.0); }
933 		if(difc!=0) difc=pow(difc,-1.0/3.0); }
934 
935 	return difc; }
936 
937 
938 /* bngmakedisplaysize */
bngmakedisplaysize(bngptr bng,int index,int totalmn)939 double bngmakedisplaysize(bngptr bng,int index,int totalmn) {
940 	double displaysize;
941 	int mn,j;
942 	simptr sim;
943 	enum MolecState ms;
944 
945 	displaysize=0;
946 	sim=bng->bngss->sim;
947 	j=stringfind(sim->mols->spname,sim->mols->nspecies,bng->bspshortnames[index]);	// look for species in Smoldyn species list
948 
949 	if(j>0) {
950 		ms=bng->bspstate[index];
951 		if(ms==MSbsoln) ms=MSsoln;
952 		displaysize=sim->mols->display[j][ms]; }
953 	else if(totalmn==1) {      // just 1 monomer
954 		for(mn=0;mn<bng->nmonomer;mn++)
955 			if(bng->monomercount[mn]>0) {
956 				displaysize=bng->monomerdisplaysize[mn];
957 				mn=bng->nmonomer; }}
958 	else {                // multiple monomers
959 		for(mn=0;mn<bng->nmonomer;mn++)
960 			if(bng->monomercount[mn]>0) {
961 				displaysize+=bng->monomercount[mn]*pow(bng->monomerdisplaysize[mn],3.0); }
962 		if(displaysize>0) displaysize=pow(displaysize,1.0/3.0); }
963 
964 	return displaysize; }
965 
966 
967 /* bngmakecolor */
bngmakecolor(bngptr bng,int index,int totalmn,double * color)968 int bngmakecolor(bngptr bng,int index,int totalmn,double *color) {
969 	double weight,totalweight;
970 	int mn,j;
971 	simptr sim;
972 	enum MolecState ms;
973 
974 	color[0]=color[1]=color[2]=0;
975 	sim=bng->bngss->sim;
976 	j=stringfind(sim->mols->spname,sim->mols->nspecies,bng->bspshortnames[index]);	// look for species in Smoldyn species list
977 
978 	if(j>0) {
979 		ms=bng->bspstate[index];
980 		if(ms==MSbsoln) ms=MSsoln;
981 		color[0]=sim->mols->color[j][ms][0];
982 		color[1]=sim->mols->color[j][ms][1];
983 		color[2]=sim->mols->color[j][ms][2]; }
984 	else if(totalmn==1) {      // just 1 monomer
985 		for(mn=0;mn<bng->nmonomer;mn++)
986 			if(bng->monomercount[mn]>0) {
987 				color[0]=bng->monomercolor[mn][0];
988 				color[1]=bng->monomercolor[mn][1];
989 				color[2]=bng->monomercolor[mn][2];
990 				mn=bng->nmonomer; }}
991 	else {                // multiple monomers
992 		weight=0;
993 		totalweight=0;
994 		for(mn=0;mn<bng->nmonomer;mn++)
995 			if(bng->monomercount[mn]>0) {
996 				weight=bng->monomercount[mn]*bng->monomerdisplaysize[mn];
997 				totalweight+=weight;
998 				color[0]+=weight*bng->monomercolor[mn][0];
999 				color[1]+=weight*bng->monomercolor[mn][1];
1000 				color[2]+=weight*bng->monomercolor[mn][2]; }
1001 		color[0]/=totalweight;
1002 		color[1]/=totalweight;
1003 		color[2]/=totalweight; }
1004 	return 0; }
1005 
1006 
1007 /* bngmakesurfaction */
bngmakesurfaction(bngptr bng,int index,int totalmn,enum SrfAction ** srfaction,surfactionptr ** actdetails)1008 void bngmakesurfaction(bngptr bng,int index,int totalmn,enum SrfAction **srfaction,surfactionptr **actdetails) {
1009 	int j,s,mn;
1010 	simptr sim;
1011 	enum MolecState ms;
1012 	enum SrfAction trialact,currentact;
1013 	enum PanelFace face;
1014 	surfaceptr srf;
1015 	surfactionptr trialdet,currentdet;
1016 
1017 	sim=bng->bngss->sim;
1018 	j=stringfind(sim->mols->spname,sim->mols->nspecies,bng->bspshortnames[index]);	// look for species in Smoldyn species list
1019 	for(s=0;s<bng->bngmaxsurface;s++)
1020 		srfaction[s][PFfront]=srfaction[s][PFback]=SAtrans;
1021 
1022 	ms=bng->bspstate[index];
1023 	if(ms==MSbsoln) ms=MSsoln;
1024 
1025 	if(j>0) {
1026 		for(s=0;s<bng->bngmaxsurface;s++) {
1027 			srf=sim->srfss->srflist[s];
1028 			for(face=PFfront;face<=PFnone;face=(enum PanelFace)(face+1)) {
1029 				srfaction[s][face]=srf->action[j][ms][face];
1030 				actdetails[s][face]=srf->actdetails[j][ms][face]; }}}
1031 
1032 	else if(totalmn==1) {
1033 		for(mn=0;mn<bng->nmonomer;mn++) {
1034 			if(bng->monomercount[mn]>0) {
1035 				for(s=0;s<bng->bngmaxsurface;s++) {
1036 					for(face=PFfront;face<=PFnone;face=(enum PanelFace)(face+1)) {
1037 						srfaction[s][face]=bng->monomeraction[mn][s][face];
1038 						actdetails[s][face]=bng->monomeractdetails[mn][s][face]; }}
1039 				mn=bng->nmonomer; }}}
1040 
1041 	else {
1042 		for(s=0;s<bng->bngmaxsurface;s++)
1043 			for(face=PFfront;face<=PFnone;face=(enum PanelFace)(face+1))
1044 				for(mn=0;mn<bng->nmonomer;mn++)
1045 					if(bng->monomercount[mn]>0) {	// priority: SAtrans < SAmult < SAreflect < SAjump < SAabsorb < SAport
1046 						trialact=bng->monomeraction[mn][s][face];
1047 						trialdet=bng->monomeractdetails[mn][s][face];
1048 						currentact=srfaction[s][face];
1049 						currentdet=actdetails[s][face];
1050 						if(srfcompareaction(currentact,currentdet,trialact,trialdet)>0) {
1051 							srfaction[s][face]=trialact;
1052 							actdetails[s][face]=trialdet; }}}
1053 
1054 	return; }
1055 
1056 
1057 /* bngparsespecies */
bngparsespecies(bngptr bng,int index)1058 int bngparsespecies(bngptr bng,int index) {
1059   int i1,i2,i3,mn,totalmn,er,hasmods,s;
1060   char *longname,ch,ch3,monomername[STRCHAR];
1061   simptr sim;
1062   double value,poslo[DIMMAX],poshi[DIMMAX],difc,displaysize,color[3];
1063 	enum MolecState ms,ms0,ms1,ms2,ms4;
1064 	enum SrfAction **action;
1065 	surfactionptr **actdetails;
1066 	surfaceptr srf;
1067 	enum PanelFace face;
1068 
1069   sim=bng->bngss->sim;
1070 
1071   for(mn=0;mn<bng->nmonomer;mn++)           // clear monomer counts
1072     bng->monomercount[mn]=0;
1073 
1074   longname=bng->bsplongnames[index];        // count numbers of each monomer in the longname
1075   totalmn=0;
1076 	hasmods=0;
1077   i1=i2=i3=0;
1078   while(longname[i2]!='\0') {			// i1 and i2 walk along string, i2 is ahead, they bracket monomer names
1079     while((ch=longname[i2])!='(' && ch!='.' && ch!='\0') i2++;
1080 
1081 		if(ch=='(') {									// i3 and ch3 look for a ~ character to test for modifications
1082 			i3=i2;
1083 			while((ch3=longname[i3])!='\0' && ch3!=')' && ch3!='~') i3++;
1084 			if(ch3=='~') hasmods=1; }
1085 
1086 		strncpy(monomername,longname+i1,i2-i1);
1087 		monomername[i2-i1]='\0';
1088     mn=stringfind(bng->monomernames,bng->nmonomer,monomername);
1089     if(mn<0) {
1090       mn=bngaddmonomer(bng,monomername,MSsoln);
1091       if(mn<0) return -1; }
1092     bng->monomercount[mn]++;
1093     totalmn++;
1094     if(ch=='.') i1=i2=i2+1;
1095     else if(ch=='(') {
1096       i1=i2=strparenmatch(longname,i2)+1;
1097       if(i1<=0) return -2;		// cannot parse name
1098 			while(longname[i1]=='.') i1=i2=i1+1; }}
1099 
1100 	bngmakeshortname(bng,index,totalmn,hasmods);
1101 	ms=bngmakedefaultstate(bng,index,totalmn);
1102 	difc=bngmakedifc(bng,index,totalmn);
1103 	displaysize=bngmakedisplaysize(bng,index,totalmn);
1104 	bngmakecolor(bng,index,totalmn,color);
1105 	action=NULL;
1106 	actdetails=NULL;
1107 	if(bng->bngmaxsurface) {
1108 		action=(enum SrfAction**) calloc(bng->bngmaxsurface,sizeof(enum SrfAction*));
1109 		actdetails=(surfactionptr **) calloc(bng->bngmaxsurface,sizeof(surfactionptr*));
1110 		if(!action || !actdetails) return -1;
1111 		for(s=0;s<bng->bngmaxsurface;s++) {
1112 			action[s]=(enum SrfAction*) calloc(3,sizeof(enum SrfAction));
1113 			actdetails[s]=(surfactionptr*) calloc(3,sizeof(surfactionptr));
1114 			if(!action[s] || !actdetails[s]) return -1;
1115 			action[s][PFfront]=action[s][PFback]=action[s][PFnone]=SAtrans;
1116 			actdetails[s][PFfront]=actdetails[s][PFback]=actdetails[s][PFnone]=NULL; }
1117 		bngmakesurfaction(bng,index,totalmn,action,actdetails); }
1118 
1119   i1=moladdspecies(sim,bng->bspshortnames[index]);            // add to Smoldyn simulation
1120   if(i1==-1) return -1;         // out of memory
1121   else if(i1==-4) return -3;    // illegal name
1122   else if(i1==-5)               // already exists
1123     i1=stringfind(sim->mols->spname,sim->mols->nspecies,bng->bspshortnames[index]);
1124   else if(i1==-6) return -3;    // illegal name
1125 
1126   bng->spindex[index]=i1;												// set spindex element
1127 
1128 	if(ms==MSbsoln) ms=MSsoln;
1129 	molsetdifc(sim,i1,NULL,ms,difc);							// set diffusion coefficient
1130 	molsetdisplaysize(sim,i1,NULL,MSall,displaysize);
1131 	molsetcolor(sim,i1,NULL,MSall,color);
1132 	if(bng->bngmaxsurface) {
1133 		for(s=0;s<bng->bngmaxsurface;s++) {
1134 			srf=sim->srfss->srflist[s];
1135 			for(face=(enum PanelFace)0;face<=PFnone;face=(enum PanelFace)(face+1)) {
1136 				surfsetaction(srf,i1,NULL,ms,face,action[s][face],-1);
1137 				if(action[s][face]==SAmult) {
1138 					for(ms4=(enum MolecState)0;ms4<(enum MolecState)MSMAX1;ms4=(enum MolecState)(ms4+1)) {
1139 						srfindex2tristate(ms,face,ms4,&ms0,&ms1,&ms2);
1140 						surfsetrate(srf,i1,NULL,ms0,ms1,ms2,i1,actdetails[s][face]->srfrate[ms4],1); }}}}
1141 		for(s=0;s<bng->bngmaxsurface;s++) {
1142 			free(action[s]);
1143 			free(actdetails[s]); }
1144 		free(action);
1145 		free(actdetails); }
1146 
1147   if(bng->bspcountstr[index]) {                 // parse count information
1148     value=strmatheval(bng->bspcountstr[index],bng->paramnames,bng->paramvalues,bng->nparams);
1149     if(strmatherror(NULL,1)) return -4;   // cannot parse count string
1150     bng->bspcount[index]=value;
1151     if(value>0.5) {                             // add molecule to Smoldyn simulation
1152       systemcorners(sim,poslo,poshi);
1153 			if(ms==MSsoln)
1154 				er=addmol(sim,(int)(value+0.5),bng->spindex[index],poslo,poshi,0);
1155 			else {
1156 				er=addsurfmol(sim,(int)(value+0.5),bng->spindex[index],bng->bspstate[index],NULL,NULL,-1,PSall,NULL); }
1157       if(er==1) return -1;
1158 			else if(er==2) return -6;			// no surface panels
1159       else if(er==3) return -5; }}  // too many molecules
1160 
1161   return 0; }
1162 
1163 
1164 /* bngaddspecies */
bngaddspecies(bngptr bng,int bindex,const char * longname,const char * countstr)1165 int bngaddspecies(bngptr bng,int bindex,const char *longname,const char *countstr) {
1166   int er;
1167 
1168 	if(bindex>=bng->maxbspecies) {		// enlarge list if needed
1169 		bng=bngalloc(bng,0,0,2*bindex+1,0);
1170 		if(!bng) return -1; }
1171 
1172 	if(longname) strncpy(bng->bsplongnames[bindex],longname,STRCHAR-1);
1173   else bng->bsplongnames[bindex][0]='\0';
1174 	if(countstr) strncpy(bng->bspcountstr[bindex],countstr,STRCHAR-1);
1175   else bng->bspcountstr[bindex][0]='\0';
1176 	if(bng->nbspecies<=bindex)
1177 		bng->nbspecies=bindex+1;
1178   er=bngparsespecies(bng,bindex);
1179   return er; }
1180 
1181 
1182 /******************************************************************************/
1183 /************************* structure set up - reactions ***********************/
1184 /******************************************************************************/
1185 
1186 
1187 /* bngparsereaction */
bngparsereaction(bngptr bng,int index)1188 int bngparsereaction(bngptr bng,int index) {
1189   int i1,i2,order,nprod,react[2],prod[2],er,allsoln;
1190   char string[STRCHAR];
1191   enum MolecState rctstate[2],prdstate[2];
1192 	double rate;
1193   simptr sim;
1194   rxnptr rxn;
1195 
1196   sim=bng->bngss->sim;
1197   order=sscanf(bng->brxnreactstr[index],"%i,%i",&i1,&i2);   // reactant list
1198   bng->brxnorder[index]=order;
1199   if(order>=1) {
1200     bng->brxnreact[index][0]=i1;
1201     react[0]=bng->spindex[i1];
1202 		rctstate[0]=bng->bspstate[i1]; }
1203   else {
1204     bng->brxnreact[index][0]=0;
1205     react[0]=0;
1206 		rctstate[0]=MSsoln; }
1207   if(order==2) {
1208     bng->brxnreact[index][1]=i2;
1209     react[1]=bng->spindex[i2];
1210 		rctstate[1]=bng->bspstate[i2]; }
1211   else {
1212     bng->brxnreact[index][1]=0;
1213     react[1]=0;
1214 		rctstate[1]=MSsoln; }
1215 
1216   nprod=sscanf(bng->brxnprodstr[index],"%i,%i",&i1,&i2);    // product list
1217   bng->brxnnprod[index]=nprod;
1218   if(nprod>=1) {
1219     bng->brxnprod[index][0]=i1;
1220     prod[0]=bng->spindex[i1];
1221 		prdstate[0]=bng->bspstate[i1]; }
1222   else {
1223     bng->brxnprod[index][0]=0;
1224     prod[0]=0;
1225 		prdstate[0]=MSsoln; }
1226   if(nprod==2) {
1227     bng->brxnprod[index][1]=i2;
1228     prod[1]=bng->spindex[i2];
1229 		prdstate[1]=bng->bspstate[i2]; }
1230   else {
1231     bng->brxnprod[index][1]=0;
1232     prod[1]=0;
1233 		prdstate[1]=MSsoln; }
1234 
1235 	allsoln=1;
1236 	if(order>=1 && (rctstate[0]==MSup || rctstate[0]==MSdown || rctstate[0]==MSfront || rctstate[0]==MSback)) allsoln=0;
1237 	if(order==2 && (rctstate[1]==MSup || rctstate[1]==MSdown || rctstate[1]==MSfront || rctstate[1]==MSback)) allsoln=0;
1238 	if(allsoln) {
1239 		rctstate[0]=rctstate[1]=MSsoln;
1240 		prdstate[0]=prdstate[1]=MSsoln; }
1241 
1242   snprintf(string,STRCHAR,"%s_%i",bng->bngname,index);     // reaction name
1243 
1244   rxn=RxnAddReaction(sim,string,order,react,rctstate,nprod,prod,prdstate,NULL,NULL);
1245   if(!rxn) return 1;
1246   bng->brxn[index]=rxn;
1247 
1248 	rate=strmatheval(bng->brxnratestr[index],bng->paramnames,bng->paramvalues,bng->nparams);
1249 	if(strmatherror(NULL,1)) return 2;   // cannot parse rate string
1250 	rate*=(order==1)?bng->unirate:bng->birate;
1251 	if(order==2 && react[0]==react[1]) rate*=2;	// BioNetGen divides rate by 2 so undo it here
1252 	er=RxnSetValue(sim,"rate",rxn,rate);
1253 	if(er==4) return 2;
1254   return 0; }
1255 
1256 
1257 /* bngaddreaction */
bngaddreaction(bngptr bng,int bindex,const char * reactants,const char * products,const char * rate)1258 int bngaddreaction(bngptr bng,int bindex,const char *reactants,const char *products,const char *rate) {
1259   int er;
1260 
1261 	if(bindex>=bng->maxbrxns) {		// enlarge list if needed
1262 		bng=bngalloc(bng,0,0,0,2*bindex+1);
1263 		if(!bng) return 1; }
1264 
1265 	if(reactants) strcpy(bng->brxnreactstr[bindex],reactants);
1266   else bng->brxnreactstr[bindex][0]='\0';
1267 	if(products) strcpy(bng->brxnprodstr[bindex],products);
1268   else bng->brxnprodstr[bindex][0]='\0';
1269 	if(rate) strcpy(bng->brxnratestr[bindex],rate);
1270   else bng->brxnratestr[bindex][0]='\0';
1271 	if(bng->nbrxns<=bindex)
1272 		bng->nbrxns=bindex+1;
1273   er=bngparsereaction(bng,bindex);
1274   return er; }
1275 
1276 
1277 /******************************************************************************/
1278 /************************** structure set up - groups *************************/
1279 /******************************************************************************/
1280 
1281 
1282 /* bngaddgroup */
bngaddgroup(bngptr bng,int gindex,const char * gname,const char * specieslist)1283 int bngaddgroup(bngptr bng,int gindex,const char *gname,const char *specieslist) {
1284 	int er,itct,imol,i1;
1285 	const char *splistptr;
1286 
1287 	(void)gindex;
1288 	er=moladdspeciesgroup(bng->bngss->sim,gname,NULL,0);
1289 	if(er) return 1;
1290 	splistptr=specieslist;
1291 	itct=sscanf(splistptr,"%i",&i1);
1292 	while(itct) {
1293 		imol=bng->spindex[i1];
1294 		er=moladdspeciesgroup(bng->bngss->sim,gname,NULL,imol);
1295 		if(er) return 1;
1296 		splistptr=strchr(splistptr,',');
1297 		if(splistptr) {
1298 			splistptr++;
1299 			itct=sscanf(splistptr,"%i",&i1); }
1300 		else itct=0; }
1301 
1302 	return 0; }
1303 
1304 
1305 /******************************************************************************/
1306 /****************** structure set up - high level functions *******************/
1307 /******************************************************************************/
1308 
1309 
1310 /* bngrunBNGL2 */
bngrunBNGL2(bngptr bng,char * filename,char * outname)1311 int bngrunBNGL2(bngptr bng,char *filename,char *outname) {
1312 	char string[STRCHAR],*dot;
1313 	FILE *fptr;
1314 	int vflag;
1315 
1316 	vflag=strchr(bng->bngss->sim->flags,'v')?1:0;
1317 
1318 	fptr=fopen(bng->bngss->BNG2path,"r");	// check for software
1319 	if(!fptr) return 1;
1320 	fclose(fptr);
1321 
1322 	fptr=fopen(filename,"r");		// check for input file
1323 	if(!fptr) return 2;
1324 	fclose(fptr);
1325 
1326 	strcpy(outname,filename);		// create outname
1327 	dot=strrchr(outname,'.');
1328 	if(!dot) dot=outname+strlen(outname);
1329 	strcpy(dot,".net");
1330 	remove(outname);						// delete output file
1331 
1332 	snprintf(string,STRCHAR,"perl -v > %s",outname);	// test for perl
1333 	system(string);
1334 	fptr=fopen(outname,"r");
1335 	if(!fptr) return 4;
1336 	remove(outname);
1337 
1338 	snprintf(string,STRCHAR,"perl %s %s %s",bng->bngss->BNG2path,filename,vflag?"":DEVNULL);
1339 	simLog(bng->bngss->sim,2," Running BNG2.pl on %s\n",filename);
1340 	system(string);							// generate network
1341 
1342 	fptr=fopen(outname,"r");		// check for output file
1343 	if(!fptr) return 3;					// output file was not written
1344 	fclose(fptr);
1345 
1346 	simLog(bng->bngss->sim,2," BNG2.pl ran successfully\n");
1347 	return 0; }
1348 
1349 
1350 /* bngreadstring */
bngreadstring(simptr sim,ParseFilePtr pfp,bngptr bng,const char * word,char * line2)1351 bngptr bngreadstring(simptr sim,ParseFilePtr pfp,bngptr bng,const char *word,char *line2) {
1352 	char str1[STRCHAR],str2[STRCHAR],str3[STRCHAR];
1353 	bngssptr bngss;
1354 	int itct,index,er;
1355 	enum MolecState ms;
1356 	static int inparams=0,inspecies=0,inrxns=0,ingroups=0;
1357 	double f1,v2[4];
1358 
1359 	int nvar;
1360 	char **varnames;
1361 	double *varvalues;
1362 
1363 	nvar=sim->nvar;
1364 	varnames=sim->varnames;
1365 	varvalues=sim->varvalues;
1366 
1367 	bngss=sim->bngss;
1368 
1369 	if(inparams) {																// parameters
1370 		if(!strcmp(word,"end")) inparams=0;
1371 		else {
1372 			itct=sscanf(line2,"%s %s",str1,str2);
1373 			CHECKS(itct==2,"failed to read parameter line");
1374 			CHECKS(bngaddparameter(bng,str1,str2)>=0,"failed to add parameter"); }}
1375 
1376 	else if(inspecies) {													// species
1377 		if(!strcmp(word,"end")) inspecies=0;
1378 		else {
1379 			itct=sscanf(word,"%i",&index);
1380 			CHECKS(itct==1,"failed to read species index value");
1381 			itct=sscanf(line2,"%s %s",str1,str2);
1382 			CHECKS(itct>=1,"failed to read species on species line");
1383 			if(itct==1) strcpy(str2,"0");
1384 			CHECKS(bngaddspecies(bng,index,str1,str2)==0,"failed to add species"); }}
1385 
1386 	else if(inrxns) {															// reactions
1387 		if(!strcmp(word,"end")) inrxns=0;
1388 		else {
1389 			itct=sscanf(word,"%i",&index);
1390 			CHECKS(itct==1,"failed to read reaction index value");
1391 			itct=sscanf(line2,"%s %s %s",str1,str2,str3);
1392 			CHECKS(itct==3,"failed to read all 3 values of reaction line");
1393 			CHECKS(bngaddreaction(bng,index,str1,str2,str3)==0,"failed to add reaction"); }}
1394 
1395 	else if(ingroups) {														// groups
1396 		if(!strcmp(word,"end")) ingroups=0;
1397 		else {
1398 			itct=sscanf(word,"%i",&index);
1399 			CHECKS(itct==1,"failed to read group index value");
1400 			itct=sscanf(line2,"%s %s",str1,str2);
1401 			CHECKS(itct==2,"failed to read all 2 values of groups line");
1402 			CHECKS(bngaddgroup(bng,index,str1,str2)==0,"failed to add group"); }}
1403 
1404 	else if(!strncmp(word,"setOption",9));				// setOption... (intentionally ignored)
1405 
1406 	else if(!strcmp(word,"begin")) {							// begin
1407     CHECKS(bng,"need to enter a name for this set of BNG complexes before entering its data");
1408 		itct=sscanf(line2,"%s",str1);
1409 		CHECKS(itct==1,"error reading because 'begin' needs to followed by 'parameters', 'species', or 'reactions'");
1410 		if(!strcmp(str1,"parameters")) inparams=1;
1411 		else if(!strcmp(str1,"species")) inspecies=1;
1412 		else if(!strcmp(str1,"reactions")) inrxns=1;
1413 		else if(!strcmp(str1,"groups")) ingroups=1;
1414 		CHECKS(!strnword(line2,2),"unexpected text at end of begin statement"); }
1415 
1416 	else if(!strcmp(word,"name")) {								// name
1417 		itct=sscanf(line2,"%s",str1);
1418 		CHECKS(itct==1,"error reading bng name");
1419 		bng=bngaddbng(sim,str1);
1420 		CHECKS(bng,"failed to add bng structure");
1421 		CHECKS(!strnword(line2,2),"unexpected text following name"); }
1422 
1423 	else if(!strcmp(word,"BNG2_path")) {					// BNG2_path
1424 		bngsetBNG2path(bng,line2); }
1425 
1426 	else if(!strcmp(word,"multiply")) {						// multiply
1427 		itct=strmathsscanf(line2,"%s %mlg",varnames,varvalues,nvar,str1,&f1);
1428 		CHECKS(itct==2,"multiply format: parameter amount");
1429 		er=bngsetparam(bng,str1,f1);
1430 		CHECKS(er!=1,"unrecognized parameter");
1431 		CHECKS(er!=2,"multiply amount needs to be greater than 0");
1432 		CHECKS(!strnword(line2,3),"unrecognized text following multiply statement"); }
1433 
1434 	else if(!strcmp(word,"monomer") || !strcmp(word,"monomers")) {	// monomer or monomers
1435 		while((itct=sscanf(line2,"%s",str1))==1) {
1436 			er=bngaddmonomer(bng,str1,MSsoln);
1437 			CHECKS(er!=-1,"out of memory adding monomer %s",str1);
1438 			CHECKS(er!=-2,"monomer name '%s' is not permitted",str1);
1439 			CHECKS(!er,"error adding monomer %s",str1); }}
1440 
1441 	else if(!strcmp(word,"monomer_difc")) {				// monomer_difc
1442 		itct=strmathsscanf(line2,"%s %mlg",varnames,varvalues,nvar,str1,&f1);
1443 		CHECKS(itct==2,"monomer_difc format: monomer difc");
1444 		CHECKS(f1>=0,"monomer diffusion coefficient has to be >=0");
1445 		er=bngsetmonomerdifc(bng,str1,f1);
1446 		CHECKS(er!=-1,"out of memory adding monomer %s",str1);
1447 		CHECKS(er!=-2,"invalid '%s' monomer name",str1);
1448 		CHECKS(!er,"failed to set monomer diffusion coefficient");
1449 		CHECKS(!strnword(line2,3),"unexpected text following monomer_difc"); }
1450 
1451 	else if(!strcmp(word,"monomer_display_size")) {	// monomer_display_size
1452 		itct=strmathsscanf(line2,"%s %mlg",varnames,varvalues,nvar,str1,&f1);
1453 		CHECKS(itct==2,"monomer_display_size format: monomer display_size");
1454 		er=bngsetmonomerdisplaysize(bng,str1,f1);
1455 		CHECKS(er!=-1,"out of memory adding monomer %s",str1);
1456 		CHECKS(er!=-2,"invalid '%s' monomer name",str1);
1457 		CHECKS(!er,"failed to set monomer display size");
1458 		CHECKS(!strnword(line2,3),"unexpected text following monomer_display_size"); }
1459 
1460 	else if(!strcmp(word,"monomer_color") || !strcmp(word,"monomer_colour")) {				// monomer_color
1461 		itct=sscanf(line2,"%s",str1);
1462 		CHECKS(itct==1,"monomer_color format: monomer color");
1463 		line2=strnword(line2,2);
1464 		CHECKS(line2,"monomer_color format: monomer color");
1465 		er=graphicsreadcolor(&line2,v2);
1466 		CHECKS(er!=3,"color values need to be between 0 and 1");
1467 		CHECKS(er!=4,"color name not recognized");
1468 		CHECKS(er!=6,"alpha values need to be between 0 and 1");
1469 		CHECKS(er==0,"format is either 3 numbers or color name");
1470 		er=bngsetmonomercolor(bng,str1,v2);
1471 		CHECKS(er!=-1,"out of memory adding monomer %s",str1);
1472 		CHECKS(er!=-2,"invalid '%s' monomer name",str1);
1473 		CHECKS(!er,"failed to set monomer display size");
1474 		CHECKS(!strnword(line2,3),"unexpected text following monomer_display_size"); }
1475 
1476 	else if(!strcmp(word,"monomer_state")) {			// monomer_state
1477 		itct=sscanf(line2,"%s %s",str1,str2);
1478 		CHECKS(itct==2,"monomer_state format: monomer state");
1479 		ms=molstring2ms(str2);
1480 		CHECKS(ms<MSMAX1,"monomer state name not allowed or not recognized");
1481 		er=bngsetmonomerstate(bng,str1,ms);
1482 		CHECKS(er!=-1,"out of memory adding monomer %s",str1);
1483 		CHECKS(er!=-2,"invalid '%s' monomer name",str1);
1484 		CHECKS(!er,"failed to set monomer default state");
1485 		CHECKS(!strnword(line2,3),"unexpected text following monomer_state"); }
1486 
1487 	else if(!strcmp(word,"expand_rules")) {				// expand_rules
1488 		itct=sscanf(line2,"%s",str1);
1489 		CHECKS(itct==1,"error reading file name");
1490 		er=bngrunBNGL2(bng,str1,str2);
1491 		CHECKS(er!=1,"BNG2.pl software is not at '%s'",bngss->BNG2path);
1492 		CHECKS(er!=2,"BioNetGen file '%s' not found",str1);
1493 		CHECKS(er!=3,"Error reading BioNetGen file '%s'. Try again with '-v' flag for more information.",str1);
1494 		CHECKS(er!=4,"Unable to run perl. Check that it is installed.");
1495 		CHECKS(!strnword(line2,2),"unexpected text following file name"); }
1496 
1497 	else {																				// unknown word
1498 		CHECKS(0,"syntax error within bngnet block: statement not recognized"); }
1499 
1500 	return bng;
1501 
1502 failure:
1503 	simParseError(sim,pfp);
1504 	return NULL; }
1505 
1506 
1507 /* loadbng */
loadbng(simptr sim,ParseFilePtr * pfpptr,char * line2)1508 int loadbng(simptr sim,ParseFilePtr *pfpptr,char* line2) {
1509 	ParseFilePtr pfp;
1510 	char word[STRCHAR],errstring[STRCHAR];
1511 	int done,pfpcode,firstline2;
1512 	bngptr bng;
1513 
1514 	pfp=*pfpptr;
1515 	done=0;
1516 	bng=NULL;
1517 	firstline2=line2?1:0;
1518 
1519 	while(!done) {
1520 		if(pfp->lctr==0)
1521 			simLog(sim,2," Reading file: '%s'\n",pfp->fname);
1522 		if(firstline2) {
1523 			strcpy(word,"name");
1524 			pfpcode=1;
1525 			firstline2=0; }
1526 		else
1527 			pfpcode=Parse_ReadLine(&pfp,word,&line2,errstring);
1528 		*pfpptr=pfp;
1529 		CHECKS(pfpcode!=3,"%s",errstring);
1530 
1531 		if(pfpcode==0);																// already taken care of
1532 		else if(pfpcode==2) {													// end reading
1533 			done=1; }
1534 		else if(!strcmp(word,"end_bng")) {					// end_bng
1535 			CHECKS(!line2,"unexpected text following end_bng");
1536 			return 0; }
1537 		else {
1538 			bng=bngreadstring(sim,pfp,bng,word,line2);
1539 			CHECK(bng); }}															// failed but error has already been reported
1540 
1541 	CHECKS(0,"end of file encountered before end_bng statement");	// end of file
1542 
1543 failure:																					// failure
1544 	if(ErrorType!=1) simParseError(sim,pfp);
1545 	*pfpptr=pfp=NULL;
1546 	return 1; }
1547 
1548 
1549 /* bngupdateparams */
bngupdateparams(simptr sim)1550 int bngupdateparams(simptr sim) {
1551 	(void)sim;
1552   return 0; }
1553 
1554 
1555 /* bngupdatelists */
bngupdatelists(simptr sim)1556 int bngupdatelists(simptr sim) {
1557 	(void)sim;
1558   return 0; }
1559 
1560 
1561 /* bngupdate */
bngupdate(simptr sim)1562 int bngupdate(simptr sim) {
1563 	int er;
1564 	bngssptr bngss;
1565 
1566 	bngss=sim->bngss;
1567 	if(bngss) {
1568 		if(bngss->condition<=SClists) {
1569 			er=bngupdatelists(sim);
1570 			if(er) return er;
1571 			bngsetcondition(bngss,SCparams,1); }
1572 		if(bngss->condition==SCparams) {
1573 			er=bngupdateparams(sim);
1574 			if(er) return er;
1575 			bngsetcondition(bngss,SCok,1); }}
1576   return 0; }
1577 
1578 
1579 
1580 /******************************************************************************/
1581 /*************************** core simulation functions ************************/
1582 /******************************************************************************/
1583 
1584 
1585