Skip to content

Commit

Permalink
unit debugging
Browse files Browse the repository at this point in the history
  • Loading branch information
ssandrews committed Jul 13, 2024
1 parent 9d5dd8b commit 9c04b14
Show file tree
Hide file tree
Showing 9 changed files with 107 additions and 33 deletions.
Binary file modified docs/Smoldyn/SmoldynCodeDoc.pdf
Binary file not shown.
40 changes: 39 additions & 1 deletion docs/Smoldyn/SmoldynCodeDoc.tex
Original file line number Diff line number Diff line change
Expand Up @@ -5182,7 +5182,7 @@ \subsection*{Functions}

This function starts by printing the message to a string using the variable arguments functions \ttt{va\_start}, \ttt{vsprintf}, and \ttt{va\_end}. Then, it replaces any unit placeholders with the current working units, if there are any. After that, it sends output to the logging outputs if appropriate. Then it checks the flags and message importance, and sends it to the correct output if appropriate.

As a comment on units, this expects unit placeholders to start with a pipe character and then list dimensions; parentheses are optional and the word is allowed to end with various puntuation (``, . : ; + - * / $<$ $>$ = !''). Examples: ``\ttt{|L}'', ``\ttt{|L2/T}'', ``\ttt{|(L/T)}'', ``\ttt{|/T,}''. If the file has working units, then they are substituted into the placeholders, and any parentheses or terminal comma are kept. If the file doesn't have working units, then the placeholders and any parentheses are removed, but the terminal comma is still kept.
This function can output units. It expects unit placeholders to start with a pipe character and then list dimensions. Parentheses around the complete units, which are started after the pipe character, are optional. Using parentheses means that any units are shown within parentheses but lack of units means that parentheses are not shown. The string word that contains the units is allowed to end with a single character of various puntuation (``, . : ; + - * / $<$ $>$ = ! )''), which is kept in the output; however, the word needs to have a whitespace character after that punctuation character. Examples: ``\ttt{\%g|L }'', ``\ttt{\%g|L2/T, }'', ``\ttt{text |(L/T): }'', ``\ttt{|/T, }''. If the file has working units, then they are substituted into the placeholders, and any parentheses or terminal punctuation are kept. If the file doesn't have working units, then the placeholders and any parentheses are removed, but the terminal punctuation is still kept.

\begin{longtable}[c]{clll}
\ttt{importance} & meaning & example & flags and display\\
Expand Down Expand Up @@ -6202,6 +6202,44 @@ \subsection*{Control and flow - network generation}

Within the molecule superstructure, \ttt{expand} is a flag for on-the-fly rule-based modeling. It is initialized to 0 and stays that way so long as there have never been any molecules of this species, it is increased to 1 if at least one molecule of this species has been created but it has not yet been used for rule expansion, and is set to either 2 or 3 if it has been used for expansion.

% Section: Units
\section{Units}

Units were introduced in Smoldyn version 2.74. The code is surprisingly simple, and is located primarily in \ttt{strunits}, in the string2.c library. There are a few parsing tasks.

First, input occurs with the \ttt{strmathsscanf} function, which is my version of the standard library \ttt{sscanf} function. This function describes the units for the value that will be input with text that looks like \ttt{\%mlg|L}. Here, the \ttt{m} character shows that this is a special input that needs to be read by the \ttt{strmathsscanf} parser, rather than being sent on to the \ttt{sscanf} parser. The subsequent \ttt{lg} characters are the standard formatting characters for inputting a double. The vertical bar (pipe) shows that unit information follows, and the \ttt{L} shows that the unit should be a length.

The second parsing task occurs if the user enters units with a number. Continuing with the above example, the user might enter \ttt{5.3|um}. The pipe again indicates that units follow. This value again needs to be parsed, and then checked that it is consistent with what was expected.

The third parsing task occurs when Smoldyn prints out the model value in the diagnostics. Here, Smoldyn uses \ttt{simLog} which reads \ttt{\%g|L} and then outputs \ttt{5.3um}.

\subsection{Unit formatting}

Units start with a pipe character. After that, the code generally uses ``\ttt{L}'' for length, and ``\ttt{T}'' for time. Combine multiple units with either a ``\ttt{.}'' or a ``\ttt{/}'' character. If units have a power that is not equal to 1, enter it either directly after the unit name or with a carat first. Here are examples:

\begin{longtable}[c]{lll}
\ttt{L.L} & \ttt{L2} & \ttt{L\^{}2}\\
\ttt{/L} & \ttt{L-1} & \ttt{L\^{}-1}\\
\ttt{L.L/T} & \ttt{L2/T} & \ttt{L2.T-1}\\
\ttt{/L/T} & \ttt{L-1.T-1} & \ttt{L\^{}-1/T}
\end{longtable}

If a number is explicitly unitless, then use a pipe character with nothing following it. For example, the input format string ``\ttt{\%mlg|}'' is used to input a unitless value. If the user includes units, then that is an error.

For output using \ttt{simLog}, unit formatting has two additional details. First, it's allowed to enclose the units in parentheses; for example, ``\ttt{|(L2/T)}''. This will display the units in parentheses if there are units and won't display either parentheses or units if there aren't units. Second, the units format is allowed to end with exactly one puctuation character, which is kept in the output regardless of whether units are displayed or not, where the punctuation character options are: ``\ttt{, . : ; + - * / $<$ $>$ = ! )}''.


\subsection{Unit parsing control flow}

Consider inputting a diffusion coefficient using the \ttt{difc} statement. This goes through the \ttt{difc} input section of \ttt{simreadstring}, in smolsim.c, which has the line
\begin{lstlisting}
itct=strmathsscanf(line2,"%mlg|L2/T",varnames,varvalues,nvar,&flt1);
\end{lstlisting}
This line of code sends the user's input, in \ttt{line2} off to \ttt{strmathsscanf} for parsing. That function recognizes the ``\ttt{\%m}'' portion of the formatting string and then reads in the user's expression for this input as a string (variable \ttt{expression}). If the user entered a pipe character, to indicate that units follow, then \ttt{strmathsscanf} removes it and processes the two halves of the expression string separately; the first half goes off to \ttt{strmatheval} for math parsing and evaluation and the second half gets read in as a string (\ttt{dimstr}) for unit parsing. Both the ``\ttt{L2/T}'' portion of the original format string and any units entered by the user get sent off to \ttt{strunits}, with option ``\ttt{convert}'', for unit parsing.

In \ttt{strunits}, the \ttt{dimstring}, which is ``\ttt{L2/T}'' gets parsed first. To do so, it is send back into \ttt{strunits} but now with the ``\ttt{parse}'' option. That option clears the static \ttt{dimUnit} array, and then reads the units to find that the power is 2 for the first unit type (length) and -1 for the second unit type (time). It stores these in \ttt{dimUnit} so that the array is equal to \ttt{\{2,-1\}}. It also computes the conversion factor between these units and the working units and multiplies that by the user's value. Returning to the \ttt{convert} level of \ttt{strunits}, the function copies over the \ttt{dimUnit} array to \ttt{tempunit} so it won't be lost and then parses the user's units, if there are any, in the same manner. Next, it checks unit compatibility and performs conversion. If the user didn't enter units, then the function needs to know what units are being assumed for this input. To do that, it looks on the \ttt{unitStack} to get the latest unit set.



% Chapter: Smoldyn modifications
\chapter{Smoldyn modifications}
Expand Down
1 change: 1 addition & 0 deletions examples/template.txt
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@
# The model was published in Andrews (2012) Methods for Molecular Biology, 804:519.
# It executes a Michaelis-Menten reaction within and on the surface of a 2D circle.

units um ms

# Model parameters
define K_FWD 0.001 # substrate-enzyme association reaction rate
Expand Down
2 changes: 1 addition & 1 deletion source/Smoldyn/smolcomparts.c
Original file line number Diff line number Diff line change
Expand Up @@ -501,7 +501,7 @@ void compartoutput(simptr sim) {
simLog(sim,2," %i: (",k);
for(d=0;d<dim-1;d++)
simLog(sim,2,"%g|L,",cmpt->points[k][d]);
simLog(sim,2,"%g|L )\n",cmpt->points[k][d]); }
simLog(sim,2,"%g|L)\n",cmpt->points[k][d]); }
simLog(sim,2," %i logically combined compartments\n",cmpt->ncmptl);
for(cl=0;cl<cmpt->ncmptl;cl++)
simLog(sim,2," %s %s\n",compartcl2string(cmpt->clsym[cl],string),cmpt->cmptl[cl]->cname);
Expand Down
2 changes: 1 addition & 1 deletion source/Smoldyn/smolmolec.c
Original file line number Diff line number Diff line change
Expand Up @@ -2180,7 +2180,7 @@ void molssoutput(simptr sim) {

simLog(sim,2," Overall spatial resolution:");
if(maxstep==-1 || mols->condition<SCok) simLog(sim,2," not computed\n");
else simLog(sim,2," %g\n",maxstep);
else simLog(sim,2," %g|L\n",maxstep);
simLog(sim,2,"\n");
return; }

Expand Down
28 changes: 23 additions & 5 deletions source/Smoldyn/smolreact.c
Original file line number Diff line number Diff line change
Expand Up @@ -834,8 +834,26 @@ void rxnoutput(simptr sim,int order) {
if(rxn->srf) simLog(sim,2," surface: %s\n",rxn->srf->sname);
if(rxn->multiplicity>1) simLog(sim,2," elementary rate and multiplicity: %g, %i\n",rxn->rate,rxn->multiplicity);
actualrate=rxncalcrate(sim,order,r,&pgem); // actual rate constant and pgem values for this reaction
if(rxn->rate>=0) simLog(sim,2," requested and actual rate constants: %g, %g\n",rxn->rate*rxn->multiplicity,actualrate);
else simLog(sim,2," actual rate constant: %g\n",actualrate);
if(rxn->rate>=0) {
simLog(sim,2," requested and actual rate constants");
if(order==0 && dim==1) simLog(sim,2," |(/L/T):");
else if(order==0 && dim==2) simLog(sim,2," |(/L2/T):");
else if(order==0 && dim==3) simLog(sim,2," |(/L3/T):");
else if(order==1) simLog(sim,2," |(/T):");
else if(order==2 && dim==1) simLog(sim,2," |(L/T):");
else if(order==2 && dim==2) simLog(sim,2," |(L2/T):");
else if(order==2 && dim==3) simLog(sim,2," |(L3/T):");
simLog(sim,2," %g, %g\n",rxn->rate*rxn->multiplicity,actualrate); }
else {
simLog(sim,2," actual rate constant");
if(order==0 && dim==1) simLog(sim,2," |(L-1/T):");
else if(order==0 && dim==2) simLog(sim,2," |(L-2/T):");
else if(order==0 && dim==3) simLog(sim,2," |(L-3/T):");
else if(order==1) simLog(sim,2," |(/T):");
else if(order==2 && dim==1) simLog(sim,2," |(L/T):");
else if(order==2 && dim==2) simLog(sim,2," |(L2/T):");
else if(order==2 && dim==3) simLog(sim,2," |(L3/T):");
simLog(sim,2," %g\n",actualrate); }
if(pgem>=0) simLog(sim,2," geminate recombination probability: %g\n",pgem);
if(rxn->rparamt==RPconfspread) simLog(sim,2," conformational spread reaction\n");
if(rxn->tau>=0) simLog(sim,2," characteristic time: %g|T\n",rxn->tau);
Expand Down Expand Up @@ -890,10 +908,10 @@ void rxnoutput(simptr sim,int order) {
if(step>0) simLog(sim,1," step length / binding radius: %g (%s %s steps)\n",ratio,ratio>0.1 && ratio<10?"somewhat":"very",ratio>1?"large":"small");
if(step>0 && (!rev || rparamt!=RPconfspread || rparamt!=RPbounce)) { // a normal bimolecular reaction
simLog(sim,1," activation-limited reaction rate: %g\n",1/(1/rate3-1/smolmodelrate));
simLog(sim,1," dynamics are like two-radius Smoluchowski model with radius: %g\n",bindrad*chi);
simLog(sim,1," dynamics are like C&K model with gamma' and gamma: %g, %g\n",(1-chi)/chi,bindrad*(1-chi)/chi);
simLog(sim,1," dynamics are like two-radius Smoluchowski model with radius: %g|L\n",bindrad*chi);
simLog(sim,1," dynamics are like C&K model with gamma' and gamma: %g, %g|L\n",(1-chi)/chi,bindrad*(1-chi)/chi);
lambdap=(3.18243*chi-3.40864*chi*chi+0.984385*chi*chi*chi)/pow(1-chi,2.08761);
simLog(sim,1," dynamics are like Doi model with lambda' and lambda: %g, %g\n",lambdap,lambdap/bindrad*lambdap/bindrad*dsum);
simLog(sim,1," dynamics are like Doi model with lambda' and lambda: %g, %g|/T\n",lambdap,lambdap/bindrad*lambdap/bindrad*dsum);
}}

if(rxn->rparamt!=RPbounce)
Expand Down
7 changes: 4 additions & 3 deletions source/Smoldyn/smolsim.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -138,14 +138,15 @@ void simLog(simptr sim,int importance,const char* format, ...) {
while((strptr=strchr(message,'|'))) { // process unit output
strwordcpy(word,strptr,1);
wordlen=strlen(word);
if(strchr(",.:;+-*/<>=!",word[wordlen-1])) {
if(strchr(",.:;+-*/<>=!)",word[wordlen-1])) {
word[wordlen-1]='\0';
wordlen--; }
hasparen=0;
if(word[1]=='(') {
hasparen=1;
word[wordlen-1]='\0';
wordlen--; }
if(word[wordlen-1]==')') {
word[wordlen-1]='\0';
wordlen--; }}
ans=strunits(NULL,hasparen?word+2:word+1,0,hasparen?word+2:word+1,"getunits");
if(ans==1)
strMidCat(message,strptr-message,strptr-message+wordlen+(hasparen?1:0),NULL,0,0);
Expand Down
Loading

0 comments on commit 9c04b14

Please sign in to comment.