simgi 0.1.1 → 0.2
raw patch · 50 files changed
+2633/−952 lines, 50 filesdep +mersenne-random-pure64dep ~base
Dependencies added: mersenne-random-pure64
Dependency ranges changed: base
Files
- .gitignore +5/−6
- ChangeLog +11/−0
- Makefile +3/−4
- Models/brusselator.sgl +16/−12
- Models/oregonator.sgl +18/−14
- Models/volterra.sgl +13/−10
- README +1/−0
- doc/Makefile +2/−1
- doc/simgi.html +172/−107
- doc/simgi.pdf binary
- doc/simgi.rst +150/−74
- simgi.cabal +4/−4
- src/CommandLine.hs +58/−23
- src/Engine.hs +237/−80
- src/ExtraFunctions.hs +15/−6
- src/GenericModel.hs +203/−44
- src/InputCheck.hs +70/−32
- src/InputParser.hs +266/−63
- src/Messages.hs +40/−5
- src/RpnCalc.hs +22/−9
- src/RpnData.hs +51/−2
- src/RpnParser.hs +27/−17
- src/TokenParser.hs +46/−4
- src/simgi.hs +12/−10
- test/Makefile +31/−7
- test/check_event_parser/EventParserTest.hs +248/−0
- test/check_event_parser/Makefile +15/−0
- test/check_irreversible/Makefile +19/−0
- test/check_irreversible/irreversible.1.sgl +22/−0
- test/check_irreversible/irreversible.2.sgl +22/−0
- test/check_irreversible/run_test.1.sh +60/−0
- test/check_irreversible/run_test.2.sh +60/−0
- test/check_reaction_parser/Makefile +15/−0
- test/check_reaction_parser/ReactionParserTest.hs +200/−0
- test/check_reversible/Makefile +19/−0
- test/check_reversible/average.hs +38/−0
- test/check_reversible/check_deviation.hs +46/−0
- test/check_reversible/reversible.sgl +29/−0
- test/check_reversible/run_test.sh +53/−0
- test/check_rpnstack/Makefile +15/−0
- test/check_rpnstack/RpnStackTest.hs +182/−0
- test/reversible_test/Makefile +0/−19
- test/reversible_test/average.hs +0/−34
- test/reversible_test/check_deviation.hs +0/−45
- test/reversible_test/reversible.sgl +0/−25
- test/reversible_test/run_test.sh +0/−53
- test/rpnstack_test/Makefile +0/−15
- test/rpnstack_test/RpnStackTest.hs +0/−227
- test/test_helpers/TestHelpers.hs +110/−0
- test/test_helpers/make_common +7/−0
.gitignore view
@@ -1,11 +1,10 @@ syntax: glob-src/*.o-src/*.hi src/simgi-src/*.swp-test/*.o-test/*.hi-test/*.swp+*.o+*.hi+*.swp average RpnStackTest+EventParserTest+ReactionParserTest check_deviation
ChangeLog view
@@ -1,3 +1,14 @@+2009-12-01 Markus Dittrich <haskelladdict@users.sourceforge.net>++ * 0.2 released.+ * significant overhaul of SGL language spec. This will break+ backward compatibility with simgi 0.1.+ * replaced simple RNG with a 64 bit Mersenne Twister implementation.+ * added expression statements and the ability to describe+ simulation events via the "events" block.+ * added ability to specify output format.++ 2009-06-02 Markus Dittrich <haskelladdict@users.sourceforge.net> * 0.1.1 release (identical to 0.1 with only a few minor
Makefile view
@@ -1,7 +1,7 @@ # Copyright 2009 Markus Dittrich <haskelladdict@users.sourceforge.net> # Distributed under the terms of the GNU General Public License v3 -VERSION=0.1.1+VERSION=0.2 DESTDIR= prefix=/usr mandir=$(DESTDIR)$(prefix)/share/man/man1@@ -17,7 +17,7 @@ src/PrettyPrint.hs src/RpnParser.hs src/RpnCalc.hs \ src/RpnData.hs src/TokenParser.hs -all: simgi+all: simgi simgi: $(OBJECTS) ghc -i./Models -i./src $(GHC_FLAGS_RELEASE) --make src/simgi.hs@@ -33,14 +33,13 @@ make -C test - install: simgi install -d $(docdir) install -d $(bindir) install -d $(htmldir) install -m 0755 src/simgi $(bindir)/ install -m 0644 ChangeLog COPYING AUTHORS $(docdir)/- install -m 0644 doc/*.pdf doc/*.rst $(docdir)/+ install -m 0644 doc/*.pdf $(docdir)/ install -m 0644 doc/*.html $(htmldir)/
Models/brusselator.sgl view
@@ -6,21 +6,25 @@ ------------------------------------------------------} def parameters- time = 50.0- outputIter = 50000- outputFreq = 200- systemVol = nil -- interpret rates as propensities- outputFile = "brusselator_output.dat"+ time = 50.0+ outputBuffer = 50000+ outputFreq = 200+ systemVol = nil -- interpret rates as propensities+ outputFile = "brusselator_output.dat" end def molecules- x 1000- y 2000+ x = 1000+ y = 2000 end def reactions- nil -> x { 5000 }- x -> y { 50.0 }- 2x + y -> 3x { 0.00005 }- x -> nil { 5.0 }-end + nil -> x | 5000 |+ x -> y | 50.0 |+ 2x + y -> 3x | 0.00005 |+ x -> nil | 5.0 |+end++def output+ [x,y]+end
Models/oregonator.sgl view
@@ -6,23 +6,27 @@ ------------------------------------------------------} def parameters- time = 50.0- outputIter = 50000- outputFreq = 500- systemVol = nil -- interpret rates as propensities- outputFile = "oregonator_output.dat"+ time = 50.0+ outputBuffer = 50000+ outputFreq = 500+ systemVol = nil -- interpret rates as propensities+ outputFile = "oregonator_output.dat" end def molecules- x 500- y 1000- z 2000 + x = 500+ y = 1000+ z = 2000 end def reactions- y -> x { 2.0 }- x + y -> nil { 0.1 }- x -> 2x + z { 104.0 }- 2x -> nil { 0.016 }- z -> y { 26.0 }-end + y -> x | 2.0 |+ x + y -> nil | 0.1 |+ x -> 2x + z | 104.0 |+ 2x -> nil | 0.016 |+ z -> y | 26.0 |+end++def output+ [x,y,z]+end
Models/volterra.sgl view
@@ -6,21 +6,24 @@ ------------------------------------------------------} def parameters- time = 50.0- outputIter = 50000- outputFreq = 200- systemVol = nil -- interpret rates as propensities- outputFile = "volterra_output.dat"+ time = 50.0+ outputBuffer = 50000+ outputFreq = 200+ systemVol = nil -- interpret rates as propensities+ outputFile = "volterra_output.dat" end def molecules- x 1000- y 2000+ x = 1000+ y = 2000 end def reactions- x -> 2x { 10.0 }- x + y -> 2y { 0.01 }- y -> nil { 10.0 }+ x -> 2x | 10.0 |+ x + y -> 2y | 0.01 |+ y -> nil | 10.0 | end +def output+ [x,y]+end
README view
@@ -1,6 +1,7 @@ DESCRIPTION: ------------ + simgi is a stochastic simulator using the Gillespie algorithm [1,2]. simgi is released under the GPL v3.
doc/Makefile view
@@ -21,11 +21,12 @@ clean_old:- rm -f *.html *.pdf+ rm -f *.html *.pdf .PHONY: clean clean:+ rm -f *.aux *.log *.html *.pdf *.tex *.out
doc/simgi.html view
@@ -3,14 +3,14 @@ <html xmlns="http://www.w3.org/1999/xhtml" xml:lang="en" lang="en"> <head> <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />-<meta name="generator" content="Docutils 0.5: http://docutils.sourceforge.net/" />+<meta name="generator" content="Docutils 0.6: http://docutils.sourceforge.net/" /> <title>simgi - A Stochastic Gillespie Simulator for Molecular Systems</title> <meta name="author" content="Markus Dittrich" /> <style type="text/css"> /* :Author: David Goodger (goodger@python.org)-:Id: $Id: html4css1.css 5196 2007-06-03 20:25:28Z wiemann $+:Id: $Id: html4css1.css 5951 2009-05-18 18:03:10Z milde $ :Copyright: This stylesheet has been placed in the public domain. Default cascading style sheet for the HTML output of Docutils.@@ -158,12 +158,33 @@ hr.docutils { width: 75% } -img.align-left {- clear: left }+img.align-left, .figure.align-left{+ clear: left ;+ float: left ;+ margin-right: 1em } -img.align-right {- clear: right }+img.align-right, .figure.align-right {+ clear: right ;+ float: right ;+ margin-left: 1em } +.align-left {+ text-align: left }++.align-center {+ clear: both ;+ text-align: center }++.align-right {+ text-align: right }++/* reset inner alignment in figures */+div.align-right {+ text-align: left }++/* div.align-center * { */+/* text-align: left } */+ ol.simple, ul.simple { margin-bottom: 1em } @@ -217,8 +238,7 @@ pre.address { margin-bottom: 0 ; margin-top: 0 ;- font-family: serif ;- font-size: 100% }+ font: inherit } pre.literal-block, pre.doctest-block { margin-left: 2em ;@@ -296,7 +316,7 @@ <tr class="field"><th class="docinfo-name">email:</th><td class="field-body">haskelladdict at users dot sourceforge dot net</td> </tr> <tr><th class="docinfo-name">Version:</th>-<td>0.1.1 (06/02/2009)</td></tr>+<td>0.2 (12/01/2009)</td></tr> </tbody> </table> <div class="section" id="contents">@@ -309,28 +329,22 @@ <li><a class="reference internal" href="#simgi-model-generation-language-sgl">Simgi Model Generation Language (SGL)</a></li> <li><a class="reference internal" href="#example-input-files">Example Input Files</a></li> <li><a class="reference internal" href="#bugs">Bugs</a></li>-<li><a class="reference internal" href="#references">References</a></li> </ol> </div> <div class="section" id="introduction"> <h1>Introduction</h1> <p><strong>simgi</strong> is a fairly simple and straightforward stochastic simulator-based on Gillspie's <a class="footnote-reference" href="#id7" id="id1">[1]</a> direct method. <strong>simgi</strong> is implemented in+based on Gillespie's <a class="footnote-reference" href="#id4" id="id1">[1]</a> direct method. <strong>simgi</strong> is implemented in pure Haskell, command line driven and comes with a flexible simulation description language called <a class="reference internal" href="#simgi-model-generation-language-sgl">Simgi Model Generation Language (SGL)</a>.-More information is available from the <a class="reference external" href="http://sourceforge.net/projects/simgi">project summary page</a>.</p>+<strong>simgi</strong> uses a fast 64 bit implementation of the Mersenne Twister+algorithm as random number source.</p> </div> <div class="section" id="status"> <h1>Status</h1>-<p>The 0.1 release of <strong>simgi</strong> provides a fully functional simulator-but should still be treated as an alpha version since several parts-of the code are currently not fully optimal. This is particularly-true for the random number generator which presently leverages the-StdGen instance of RandomGen <a class="footnote-reference" href="#id8" id="id2">[2]</a> and is probably not sufficient for-large simulations in terms of random number quality. Later revisions-of simgi will have a more sophisticated random number generator.-Nevertheless, for small systems (such as the examples in the-<em>Models/</em> directory) the current implementation should be sufficient.</p>+<p>The 0.2 release of <strong>simgi</strong> provides a fully functional simulator+which has been tested on several model systems some of which were+fairly large.</p> </div> <div class="section" id="download"> <h1>Download</h1>@@ -342,6 +356,7 @@ <ul class="simple"> <li><a class="reference external" href="http://haskell.org/ghc/">>=ghc-6.10</a></li> <li><a class="reference external" href="http://gmplib.org/">>=gmp-4.3</a></li>+<li><a class="reference external" href="http://hackage.haskell.org/package/mersenne-random-pure64">>=mersenne-random-pure64</a></li> </ul> <p>To compile the documentation (not required), you will also need</p> <ul class="simple">@@ -350,15 +365,16 @@ </ul> <p>Building of <strong>simgi</strong> can be done either via</p> <ul class="simple">-<li>the standard <tt class="docutils literal"><span class="pre">make,</span> <span class="pre">make</span> <span class="pre">check,</span> <span class="pre">make</span> <span class="pre">install</span></tt></li>+<li>the standard <tt class="docutils literal">make, make check, make install</tt></li> <li>or via cabal</li> </ul> </div> <div class="section" id="simgi-model-generation-language-sgl"> <h1>Simgi Model Generation Language (SGL)</h1>-<p>simgi simulations are described via <a class="reference internal" href="#simgi-model-generation-language-sgl">Simgi Model Generation Language-(SGL)</a>. The corresponding simulation files typically have an <em>.sgl</em>-extension, but this is not enforced by the <strong>simgi</strong> simulation engine.</p>+<p><strong>simgi</strong> simulations are described via <a class="reference internal" href="#simgi-model-generation-language-sgl">Simgi Model Generation Language+(SGL)</a>. The corresponding simulation input files typically have an <em>.sgl</em>+extension, but this is not enforced by the <strong>simgi</strong> simulation+engine.</p> <p>A SGL file consists of zero or more descriptor blocks of the form</p> <pre class="literal-block"> def <block name>@@ -367,120 +383,190 @@ end </pre>-<p>The (but see <a class="footnote-reference" href="#id9" id="id3">[3]</a>) formatting of the input files is very flexible. In-particular, neither newlines <a class="footnote-reference" href="#id10" id="id4">[4]</a> nor extraneous whitespace matter.-Hence, the above SDL block could also be written on a single line.+<p>The formatting of the input files is very flexible. In+particular, neither newlines <a class="footnote-reference" href="#id5" id="id2">[2]</a> nor extraneous whitespace matter.+Hence, the above SGL block could have also been written on a single line. However, it is strongly recommended to stick to a consistent and "visually simple" layout to aid in "comprehending" the underlying-model.</p>+model. Also, it is important to point out that <strong>simgi</strong>'s parser is+case sensitive.</p> <p><strong>Comments</strong> can be added to the SGL file and are parsed according to the Haskell language specs</p> <ul class="simple"> <li>simple line comments begin with a <tt class="docutils literal"><span class="pre">--</span></tt> token and treat everything-until the next newline as a comment, including valid SDL commands.-Hence, SDL blocks containing line comments need to be separated by+until the next newline as a comment, including valid SGL commands.+Hence, SGL blocks containing line comments need to be separated by newlines in order to be parsed correctly.</li>-<li>block comments begin with a <tt class="docutils literal"><span class="pre">{-</span></tt> token and end with a <tt class="docutils literal"><span class="pre">-}</span></tt>+<li>block comments begin with a <tt class="docutils literal">{-</tt> token and end with a <tt class="docutils literal"><span class="pre">-}</span></tt> token. Everything within a comment block is ignored by the parser and block comments can be nested.</li> </ul>-<p>Currently, the SDL specs define the following block types with their+<p><strong>Expression Statements</strong> are an important and useful part of SGL.+<tt class="docutils literal">Expression statements</tt> are enclosed in curly braces and can contain+any mathematical expression involving doubles, the simulation time+(via the keyword <tt class="docutils literal">TIME</tt>), as well as the values of any variable or+molecule count. The values of time, molecule counts and variables+are evaluated at run time and represent the instantaneous values at the+time at which the expression is evaluated.+<tt class="docutils literal">Expressions statements</tt> can contain any+arithmetic expression involving the standard operators "+", "-", "*", "/", "^"+(exponentiation), and the mathematical functions <tt class="docutils literal">sqrt, exp, log, log2, log10, sin,+cos, tan, asin, acos, atan, sinh, cosh, tanh, asinh, acosh, atanh, acosh, atanh,+erf, erfc, abs</tt>.</p>+<p>Internally, <tt class="docutils literal">expression statements</tt> are converted into a compute stack+in RPN format which is evaluated at run-time. Even though this+procedure is fairly efficient, there is some numerical overhead+incurred at each iteration and the use of complicated rate+expressions should therefore be avoided if possible.</p>+<p>Below is a list of all SGL blocks available for describing simulations.+Presently, the order of blocks matters and should be exactly the same+in which they are described below. Several SGL blocks are+optional and are marked as such below.</p>+<p>Currently, the SGL specs define the following block types with their respective block commands and block content:</p>-<p><strong>parameter block:</strong> <em><block name> = parameters</em></p>+<p><strong>parameter block:</strong> <tt class="docutils literal"><block name></tt> = <em>parameters</em></p> <blockquote> <p>The purpose of the parameter block is to describe the global simulation parameters. The following parameters are currently supported:</p> <dl class="docutils">-<dt><em>time = <double></em></dt>+<dt><em>time</em> = <tt class="docutils literal">Double</tt></dt> <dd>Maximum simulation time in seconds. Default is 0.0 s.</dd>-<dt><em>outputIter = <Integer></em></dt>+<dt><em>outputBuffer</em> = <tt class="docutils literal">Integer</tt></dt> <dd><p class="first">Output will be kept in memory and written to the output file and-stdout every <em>outputIter</em> iterations. Larger values should result-in faster simulations but require more system memory.+stdout every <em>outputBuffer</em> iterations. Larger values should+result in faster simulations but require more system memory. Default is to write output every 10000 iterations.</p>-<p class="last">Note: <em>outputIter</em> only affects how often output is written to-the output file, not how much is being accumulated during a+<p class="last">Note: <em>outputBuffer</em> only affects how often output is written to+the output file, not how much output is actually generated during a simulation (see outputFreq parameter).</p> </dd>-<dt><em>outputFreq = <Integer></em></dt>-<dd>Frequency with which output is generated. Default is 1000.</dd>-<dt><em>systemVol = <double></em></dt>+<dt><em>outputFreq</em> = <tt class="docutils literal">Integer</tt></dt>+<dd>Iteration frequency with which output is generated. Default is every 1000+iterations. Please note that output is written to the output file in batches of+<em>outputBuffer</em>.</dd>+<dt><em>systemVol</em> = <tt class="docutils literal">Double</tt></dt> <dd>Volume of the simulation system in liters. This is needed to properly compute the reaction rates in molar units. If rates should rather be interpreted as reaction propensities (like-in <a class="footnote-reference" href="#id7" id="id5">[1]</a>) please set <em>systemVol = nil</em>. Default is a system+in <a class="footnote-reference" href="#id4" id="id3">[1]</a>) please set <em>systemVol = nil</em>. Default is a system volume of 1.0 liter.</dd>-<dt><em>outputFile = <quoted string></em>:</dt>+<dt><em>outputFile</em> = <tt class="docutils literal">Quoted String</tt></dt> <dd>Name of the output file. This is the only required parameter in the parameter section. If not given, the simulation will terminate.</dd> </dl> </blockquote>-<p><strong>molecule block:</strong> <em><block name> = molecules</em></p>+<p><strong>variable block:</strong> <tt class="docutils literal"><block name></tt> = <em>variables</em></p> <blockquote> <p>This block consist of a list of pairs of the form</p> <pre class="literal-block">-<String> <Integer>+String = <variable expression> </pre>+<p>where <tt class="docutils literal">String</tt> is the variable name, and <tt class="docutils literal"><variable expression></tt>+is either a <tt class="docutils literal">Double</tt> or an <tt class="docutils literal">expression statement</tt> as defined above.+Variables can be used in any other <tt class="docutils literal">expression statement</tt> in the+SGL file including reaction rate definitions. Please make sure to+not define a variable in terms of itself to avoid infinite recursion.</p>+</blockquote>+<p><strong>molecule block:</strong> <tt class="docutils literal"><block name></tt> = <em>molecules</em></p>+<blockquote>+<p>This block consist of a list of pairs of the form</p>+<pre class="literal-block">+String = Integer+</pre> <p>giving the name of each molecule and the number of molecules present initially. For example, the following molecule definition-block defines molecules <tt class="docutils literal"><span class="pre">A</span></tt> and <tt class="docutils literal"><span class="pre">B</span></tt> with initial numbers of+block defines molecules <tt class="docutils literal">A</tt> and <tt class="docutils literal">B</tt> with initial numbers of 100 and 200, respectively</p> <pre class="literal-block"> def molecules- A 100- B 200+ A = 100+ B = 200 end </pre>+<p><strong>NOTE</strong>: Please do not use any of the predefined mathematical+functions or defined variables (including <tt class="docutils literal">TIME</tt>) as+molecule names since this will lead to undefined behavior.</p> </blockquote>-<p><strong>reaction block</strong>: <em><block name> = reactions</em></p>+<p><strong>reaction block</strong>: <tt class="docutils literal"><block name></tt> = <em>reactions</em></p> <blockquote> <p>This block describes the reactions between molecules defined in the molecule block. Reactions are specified via</p> <pre class="literal-block">-reactants -> product { rate expression }+<reactants> -> <products> | <rate expression> | </pre>-<p>Here, <tt class="docutils literal"><span class="pre">reactants</span></tt> and <tt class="docutils literal"><span class="pre">products</span></tt> are of the form</p>+<p>Here, <tt class="docutils literal"><reactants></tt> and <tt class="docutils literal"><products></tt> are of the form</p> <pre class="literal-block">-<Integer> <String> + <Integer> <String> + .....+Integer String + Integer String + ..... </pre>-<p>In this expression, <tt class="docutils literal"><span class="pre"><String></span></tt> is the reactant or product name-as defined in the molecule block and <tt class="docutils literal"><span class="pre"><Integer></span></tt> an optional-integer specifying the stoichiometry. If <tt class="docutils literal"><span class="pre"><Integer></span></tt> is not+<p>In this expression, <tt class="docutils literal">String</tt> is a molecule name+as defined in the molecule block and <tt class="docutils literal">Integer</tt> an optional+integer specifying the stoichiometry. If <tt class="docutils literal">Integer</tt> is not explicitly given, it is assumed to be 1.</p>-<p>The reaction rate can either be a fixed value of type <tt class="docutils literal"><span class="pre"><Double></span></tt>-or else an mathematical expression involving <tt class="docutils literal"><span class="pre"><Double></span></tt>,-molecule names, and the current simulation time. Hence, <strong>simgi</strong>-rate expressions can be arbitrary complex functions of the-instantaneous simulation time and the instantaneous numbers of any-molecule in the model. The parser will interpret any string in the-rate expression as a molecule name in a case sensitive fashion,-a mathematical operator or function (see <a class="footnote-reference" href="#id11" id="id6">[5]</a> for supported-functions), or the special variable TIME which refers to the-current simulation time. Hence, do <strong>not</strong> use any of the-mathematical keywords as a molecule name; this leads to undefined-behavior.</p>-<p>Here is an example reaction block for the two molecules <tt class="docutils literal"><span class="pre">A</span></tt> and-<tt class="docutils literal"><span class="pre">B</span></tt> defined above:</p>+<p>The <tt class="docutils literal"><rate expression></tt> can either be a fixed value of type+<tt class="docutils literal">Double</tt> or an <tt class="docutils literal">expression statement</tt> as defined above.</p>+<p>Below is an example reaction block for the two molecules <tt class="docutils literal">A</tt> and+<tt class="docutils literal">B</tt> defined above:</p> <pre class="literal-block"> define reactions- 2A + B -> A { 10.0e-5 }- B -> A { 2.0e-5 * A * exp(-0.5*TIME) }+ 2A + B -> A | 10.0e-5 |+ B -> A | {2.0e-5 * A * exp(-0.5*TIME)} | end </pre>-<p>In the first reaction, 2 <tt class="docutils literal"><span class="pre">A</span></tt> molecules react with one <tt class="docutils literal"><span class="pre">B</span></tt> to-yield another <tt class="docutils literal"><span class="pre">A</span></tt> at a rate of 10.0e-5 1/(Mol s). The second-reaction describes a decay of <tt class="docutils literal"><span class="pre">B</span></tt> back to <tt class="docutils literal"><span class="pre">A</span></tt> at a rate-that is computed based on the instantaneous number of <tt class="docutils literal"><span class="pre">A</span></tt>+<p>In the first reaction, 2 <tt class="docutils literal">A</tt> molecules react with one <tt class="docutils literal">B</tt> to+yield another <tt class="docutils literal">A</tt> at a rate of 10.0e-5. The second+reaction describes a decay of <tt class="docutils literal">B</tt> back to <tt class="docutils literal">A</tt> at a rate+that is computed based on the instantaneous number of <tt class="docutils literal">A</tt> molecules present and which decays exponentially with simulation time.</p>-<p>Internally, rate expressions are converted into a compute stack-in RPN format which is evaluated at run-time. Even though this-procedure is fairly efficient, there is some numerical overhead-incurred at each iteration and the use of complicated rate-expressions should therefore be avoided if possible.</p> </blockquote>+<p><strong>event block</strong>: <tt class="docutils literal"><block name></tt> = <em>events</em></p>+<blockquote>+<p>An event block allows one to specify events which will occur during+the simulation. Each event consists of a <tt class="docutils literal"><trigger expression></tt> and+an associated set of <tt class="docutils literal"><action expressions></tt>.+Events are specified via</p>+<pre class="literal-block">+{ <trigger expression> } => { <action expression> }+</pre>+<p>Here, <tt class="docutils literal">trigger expression</tt> is of the form</p>+<pre class="literal-block">+<trigger primitive> [ <boolean operator> <trigger primitive>]+</pre>+<p>with <tt class="docutils literal"><trigger primitive></tt> defined by</p>+<pre class="literal-block">+<expression statement> relational operator <expression statement>+</pre>+<p>Each <tt class="docutils literal"><trigger primitive></tt> contains two <tt class="docutils literal">expression statements</tt>+as defined above and a <tt class="docutils literal">relational operator</tt> which can be+any of <tt class="docutils literal">>=</tt>, <tt class="docutils literal"><=</tt>, <tt class="docutils literal">==</tt>, <tt class="docutils literal">></tt>, and <tt class="docutils literal"><</tt>. Hence, each+<tt class="docutils literal"><trigger primitive></tt> evaluates to either <tt class="docutils literal">true</tt> or <tt class="docutils literal">false</tt>.</p>+<p>Several <tt class="docutils literal"><trigger primitives></tt> can be chained together via the+<tt class="docutils literal"><boolean operators></tt> <tt class="docutils literal">&&</tt> and <tt class="docutils literal">||</tt> to yield a final boolean+value of <tt class="docutils literal">true</tt> or <tt class="docutils literal">false</tt>.</p>+<p>If the <tt class="docutils literal"><trigger expression></tt> evaluates to true during an+iteration, the associated <tt class="docutils literal"><action expressions></tt> is executed+during the same timestep.</p>+<p><tt class="docutils literal"><action expression></tt> consists of a semi-colon separated list of+assignments</p>+<pre class="literal-block">+String = <assignment expression> [; String = <assignment expression>]+</pre>+<p>where <tt class="docutils literal">String</tt> is a molecule or variable name and+<tt class="docutils literal"><expression></tt> either a <tt class="docutils literal">Double</tt> or an <tt class="docutils literal">expression statement</tt>.</p>+<p><strong>NOTE</strong>: Since molecule counts are integer values assignments+to molecule counts in <tt class="docutils literal"><action expression></tt> will be converted+to an integer value via <tt class="docutils literal">floor</tt>.</p>+</blockquote>+<p><strong>output block</strong>: <tt class="docutils literal"><block name></tt> = <em>output</em></p>+<blockquote>+<p>This block consists of a simple list of variable and molecule+names that will be streamed to the output file in the same order:</p>+<pre class="literal-block">+[ name1, name2, name3, .... ]+</pre>+</blockquote> </div> <div class="section" id="example-input-files"> <h1>Example Input Files</h1>@@ -496,37 +582,16 @@ <h1>Bugs</h1> <p>Please report all bugs and feature requests to <haskelladdict at users dot sourceforge dot net>.</p>-</div>-<div class="section" id="references">-<h1>References</h1>-<table class="docutils footnote" frame="void" id="id7" rules="none">-<colgroup><col class="label" /><col /></colgroup>-<tbody valign="top">-<tr><td class="label">[1]</td><td><em>(<a class="fn-backref" href="#id1">1</a>, <a class="fn-backref" href="#id5">2</a>)</em> Daniel T. Gillespie (1977). "Exact Stochastic Simulation of Coupled Chemical Reactions". The Journal of Physical Chemistry 81 (25): 2340-2361</td></tr>-</tbody>-</table>-<table class="docutils footnote" frame="void" id="id8" rules="none">-<colgroup><col class="label" /><col /></colgroup>-<tbody valign="top">-<tr><td class="label"><a class="fn-backref" href="#id2">[2]</a></td><td><a class="reference external" href="http://hackage.haskell.org/packages/archive/random/1.0.0.1/doc/html/System-Random#globalrng.html">http://hackage.haskell.org/packages/archive/random/1.0.0.1/doc/html/System-Random#globalrng.html</a></td></tr>-</tbody>-</table>-<table class="docutils footnote" frame="void" id="id9" rules="none">-<colgroup><col class="label" /><col /></colgroup>-<tbody valign="top">-<tr><td class="label"><a class="fn-backref" href="#id3">[3]</a></td><td>Since <strong>simgi</strong> currently is an alpha version there may be fairly drastic changes to the SDL specs in future releases until the first beta release.</td></tr>-</tbody>-</table>-<table class="docutils footnote" frame="void" id="id10" rules="none">+<table class="docutils footnote" frame="void" id="id4" rules="none"> <colgroup><col class="label" /><col /></colgroup> <tbody valign="top">-<tr><td class="label"><a class="fn-backref" href="#id4">[4]</a></td><td>An exception to this rule are line comments starting with <tt class="docutils literal"><span class="pre">--</span></tt> which ingnore everything until the next newline.</td></tr>+<tr><td class="label">[1]</td><td><em>(<a class="fn-backref" href="#id1">1</a>, <a class="fn-backref" href="#id3">2</a>)</em> Daniel T. Gillespie (1977). "Exact Stochastic Simulation of Coupled Chemical Reactions". The Journal of Physical Chemistry 81 (25): 2340-2361</td></tr> </tbody> </table>-<table class="docutils footnote" frame="void" id="id11" rules="none">+<table class="docutils footnote" frame="void" id="id5" rules="none"> <colgroup><col class="label" /><col /></colgroup> <tbody valign="top">-<tr><td class="label"><a class="fn-backref" href="#id6">[5]</a></td><td>Rate expressions can contain any arithmetic expression involving the standard operators "+", "-", "*", "/", "^" (exponentiation), and the mathematical functions <tt class="docutils literal"><span class="pre">sqrt,</span> <span class="pre">exp,</span> <span class="pre">log,</span> <span class="pre">log2,</span> <span class="pre">log10,</span> <span class="pre">sin,</span> <span class="pre">cos,</span> <span class="pre">tan,</span> <span class="pre">asin,</span> <span class="pre">acos,</span> <span class="pre">atan,</span> <span class="pre">sinh,</span> <span class="pre">cosh,</span> <span class="pre">tanh,</span> <span class="pre">asinh,</span> <span class="pre">acosh,</span> <span class="pre">atanh,</span> <span class="pre">acosh,</span> <span class="pre">atanh,</span> <span class="pre">erf,</span> <span class="pre">erfc,</span> <span class="pre">abs</span></tt>.</td></tr>+<tr><td class="label"><a class="fn-backref" href="#id2">[2]</a></td><td>An exception to this rule are line comments starting with <tt class="docutils literal"><span class="pre">--</span></tt> which ingnore everything until the next newline.</td></tr> </tbody> </table> </div>
doc/simgi.pdf view
binary file changed (92412 → 245560 bytes)
doc/simgi.rst view
@@ -6,7 +6,7 @@ :email: haskelladdict at users dot sourceforge dot net -:Version: 0.1.1 (06/02/2009)+:Version: 0.2 (12/01/2009) Contents@@ -19,30 +19,25 @@ 5) `Simgi Model Generation Language (SGL)`_ 6) `Example Input Files`_ 7) Bugs_-8) References_ + Introduction ------------ **simgi** is a fairly simple and straightforward stochastic simulator -based on Gillspie's [1]_ direct method. **simgi** is implemented in +based on Gillespie's [1]_ direct method. **simgi** is implemented in pure Haskell, command line driven and comes with a flexible simulation description language called `Simgi Model Generation Language (SGL)`_.-More information is available from the `project summary page <http://sourceforge.net/projects/simgi>`_.+**simgi** uses a fast 64 bit implementation of the Mersenne Twister+algorithm as random number source. Status ------ -The 0.1 release of **simgi** provides a fully functional simulator -but should still be treated as an alpha version since several parts -of the code are currently not fully optimal. This is particularly -true for the random number generator which presently leverages the -StdGen instance of RandomGen [2]_ and is probably not sufficient for -large simulations in terms of random number quality. Later revisions -of simgi will have a more sophisticated random number generator. -Nevertheless, for small systems (such as the examples in the -*Models/* directory) the current implementation should be sufficient. +The 0.2 release of **simgi** provides a fully functional simulator +which has been tested on several model systems some of which were+fairly large. Download@@ -58,6 +53,7 @@ - `>=ghc-6.10 <http://haskell.org/ghc/>`_ - `>=gmp-4.3 <http://gmplib.org/>`_ +- `>=mersenne-random-pure64 <http://hackage.haskell.org/package/mersenne-random-pure64>`_ To compile the documentation (not required), you will also need @@ -75,9 +71,10 @@ Simgi Model Generation Language (SGL) ------------------------------------- -simgi simulations are described via `Simgi Model Generation Language -(SGL)`_. The corresponding simulation files typically have an *.sgl* -extension, but this is not enforced by the **simgi** simulation engine. +**simgi** simulations are described via `Simgi Model Generation Language +(SGL)`_. The corresponding simulation input files typically have an *.sgl* +extension, but this is not enforced by the **simgi** simulation +engine. A SGL file consists of zero or more descriptor blocks of the form @@ -89,132 +86,218 @@ end -The (but see [3]_) formatting of the input files is very flexible. In-particular, neither newlines [4]_ nor extraneous whitespace matter. -Hence, the above SDL block could also be written on a single line. +The formatting of the input files is very flexible. In+particular, neither newlines [2]_ nor extraneous whitespace matter. +Hence, the above SGL block could have also been written on a single line. However, it is strongly recommended to stick to a consistent and "visually simple" layout to aid in "comprehending" the underlying-model.+model. Also, it is important to point out that **simgi**'s parser is +case sensitive. **Comments** can be added to the SGL file and are parsed according to the Haskell language specs - simple line comments begin with a ``--`` token and treat everything - until the next newline as a comment, including valid SDL commands. - Hence, SDL blocks containing line comments need to be separated by + until the next newline as a comment, including valid SGL commands. + Hence, SGL blocks containing line comments need to be separated by newlines in order to be parsed correctly. - block comments begin with a ``{-`` token and end with a ``-}`` token. Everything within a comment block is ignored by the parser and block comments can be nested. -Currently, the SDL specs define the following block types with their +**Expression Statements** are an important and useful part of SGL.+``Expression statements`` are enclosed in curly braces and can contain+any mathematical expression involving doubles, the simulation time +(via the keyword ``TIME``), as well as the values of any variable or +molecule count. The values of time, molecule counts and variables+are evaluated at run time and represent the instantaneous values at the+time at which the expression is evaluated.+``Expressions statements`` can contain any +arithmetic expression involving the standard operators "+", "-", "*", "/", "^" +(exponentiation), and the mathematical functions ``sqrt, exp, log, log2, log10, sin, +cos, tan, asin, acos, atan, sinh, cosh, tanh, asinh, acosh, atanh, acosh, atanh, +erf, erfc, abs``.++Internally, ``expression statements`` are converted into a compute stack+in RPN format which is evaluated at run-time. Even though this+procedure is fairly efficient, there is some numerical overhead+incurred at each iteration and the use of complicated rate +expressions should therefore be avoided if possible.+++Below is a list of all SGL blocks available for describing simulations.+Presently, the order of blocks matters and should be exactly the same+in which they are described below. Several SGL blocks are +optional and are marked as such below. ++Currently, the SGL specs define the following block types with their respective block commands and block content: -**parameter block:** *<block name> = parameters*+**parameter block:** ``<block name>`` = *parameters* The purpose of the parameter block is to describe the global simulation parameters. The following parameters are currently supported: - *time = <double>*+ *time* = ``Double`` Maximum simulation time in seconds. Default is 0.0 s. - *outputIter = <Integer>*+ *outputBuffer* = ``Integer`` Output will be kept in memory and written to the output file and - stdout every *outputIter* iterations. Larger values should result - in faster simulations but require more system memory. + stdout every *outputBuffer* iterations. Larger values should + result in faster simulations but require more system memory. Default is to write output every 10000 iterations. - Note: *outputIter* only affects how often output is written to - the output file, not how much is being accumulated during a + Note: *outputBuffer* only affects how often output is written to + the output file, not how much output is actually generated during a simulation (see outputFreq parameter). - *outputFreq = <Integer>*- Frequency with which output is generated. Default is 1000.+ *outputFreq* = ``Integer``+ Iteration frequency with which output is generated. Default is every 1000+ iterations. Please note that output is written to the output file in batches of+ *outputBuffer*. - *systemVol = <double>*+ *systemVol* = ``Double`` Volume of the simulation system in liters. This is needed to properly compute the reaction rates in molar units. If rates should rather be interpreted as reaction propensities (like in [1]_) please set *systemVol = nil*. Default is a system volume of 1.0 liter. - *outputFile = <quoted string>*:+ *outputFile* = ``Quoted String`` Name of the output file. This is the only required parameter in the parameter section. If not given, the simulation will terminate. -**molecule block:** *<block name> = molecules* +**variable block:** ``<block name>`` = *variables*+ This block consist of a list of pairs of the form :: - <String> <Integer>+ String = <variable expression> + where ``String`` is the variable name, and ``<variable expression>``+ is either a ``Double`` or an ``expression statement`` as defined above.+ Variables can be used in any other ``expression statement`` in the+ SGL file including reaction rate definitions. Please make sure to+ not define a variable in terms of itself to avoid infinite recursion.+++**molecule block:** ``<block name>`` = *molecules*++ This block consist of a list of pairs of the form ::++ String = Integer+ giving the name of each molecule and the number of molecules present initially. For example, the following molecule definition block defines molecules ``A`` and ``B`` with initial numbers of 100 and 200, respectively :: def molecules- A 100- B 200+ A = 100+ B = 200 end + **NOTE**: Please do not use any of the predefined mathematical+ functions or defined variables (including ``TIME``) as + molecule names since this will lead to undefined behavior. -**reaction block**: *<block name> = reactions* +++**reaction block**: ``<block name>`` = *reactions*+ This block describes the reactions between molecules defined in the molecule block. Reactions are specified via :: - reactants -> product { rate expression }+ <reactants> -> <products> | <rate expression> | - Here, ``reactants`` and ``products`` are of the form ::+ Here, ``<reactants>`` and ``<products>`` are of the form :: - <Integer> <String> + <Integer> <String> + .....+ Integer String + Integer String + ..... - In this expression, ``<String>`` is the reactant or product name - as defined in the molecule block and ``<Integer>`` an optional - integer specifying the stoichiometry. If ``<Integer>`` is not + In this expression, ``String`` is a molecule name + as defined in the molecule block and ``Integer`` an optional + integer specifying the stoichiometry. If ``Integer`` is not explicitly given, it is assumed to be 1. - The reaction rate can either be a fixed value of type ``<Double>`` - or else an mathematical expression involving ``<Double>``, - molecule names, and the current simulation time. Hence, **simgi** - rate expressions can be arbitrary complex functions of the - instantaneous simulation time and the instantaneous numbers of any- molecule in the model. The parser will interpret any string in the - rate expression as a molecule name in a case sensitive fashion, - a mathematical operator or function (see [5]_ for supported - functions), or the special variable TIME which refers to the - current simulation time. Hence, do **not** use any of the - mathematical keywords as a molecule name; this leads to undefined- behavior.+ The ``<rate expression>`` can either be a fixed value of type + ``Double`` or an ``expression statement`` as defined above. - Here is an example reaction block for the two molecules ``A`` and + Below is an example reaction block for the two molecules ``A`` and ``B`` defined above:: define reactions- 2A + B -> A { 10.0e-5 }- B -> A { 2.0e-5 * A * exp(-0.5*TIME) }+ 2A + B -> A | 10.0e-5 |+ B -> A | {2.0e-5 * A * exp(-0.5*TIME)} | end In the first reaction, 2 ``A`` molecules react with one ``B`` to - yield another ``A`` at a rate of 10.0e-5 1/(Mol s). The second + yield another ``A`` at a rate of 10.0e-5. The second reaction describes a decay of ``B`` back to ``A`` at a rate that is computed based on the instantaneous number of ``A`` molecules present and which decays exponentially with simulation time. - Internally, rate expressions are converted into a compute stack- in RPN format which is evaluated at run-time. Even though this- procedure is fairly efficient, there is some numerical overhead- incurred at each iteration and the use of complicated rate - expressions should therefore be avoided if possible. + +**event block**: ``<block name>`` = *events* + An event block allows one to specify events which will occur during + the simulation. Each event consists of a ``<trigger expression>`` and + an associated set of ``<action expressions>``. + Events are specified via ::++ { <trigger expression> } => { <action expression> }++ Here, ``trigger expression`` is of the form ::++ <trigger primitive> [ <boolean operator> <trigger primitive>]++ with ``<trigger primitive>`` defined by ::++ <expression statement> relational operator <expression statement>++ Each ``<trigger primitive>`` contains two ``expression statements``+ as defined above and a ``relational operator`` which can be+ any of ``>=``, ``<=``, ``==``, ``>``, and ``<``. Hence, each+ ``<trigger primitive>`` evaluates to either ``true`` or ``false``.++ Several ``<trigger primitives>`` can be chained together via the + ``<boolean operators>`` ``&&`` and ``||`` to yield a final boolean+ value of ``true`` or ``false``.++ If the ``<trigger expression>`` evaluates to true during an+ iteration, the associated ``<action expressions>`` is executed + during the same timestep.++ ``<action expression>`` consists of a semi-colon separated list of + assignments ::++ String = <assignment expression> [; String = <assignment expression>]++ + where ``String`` is a molecule or variable name and + ``<expression>`` either a ``Double`` or an ``expression statement``.++ **NOTE**: Since molecule counts are integer values assignments+ to molecule counts in ``<action expression>`` will be converted+ to an integer value via ``floor``.+++**output block**: ``<block name>`` = *output*++ This block consists of a simple list of variable and molecule+ names that will be streamed to the output file in the same order::++ [ name1, name2, name3, .... ]+++ Example Input Files ------------------- @@ -234,15 +317,8 @@ <haskelladdict at users dot sourceforge dot net>. -References------------ .. [1] Daniel T. Gillespie (1977). "Exact Stochastic Simulation of Coupled Chemical Reactions". The Journal of Physical Chemistry 81 (25): 2340-2361 -.. [2] http://hackage.haskell.org/packages/archive/random/1.0.0.1/doc/html/System-Random#globalrng.html+.. [2] An exception to this rule are line comments starting with ``--`` which ingnore everything until the next newline. -.. [3] Since **simgi** currently is an alpha version there may be fairly drastic changes to the SDL specs in future releases until the first beta release. -.. [4] An exception to this rule are line comments starting with ``--`` which ingnore everything until the next newline.--.. [5] Rate expressions can contain any arithmetic expression involving the standard operators "+", "-", "*", "/", "^" (exponentiation), and the mathematical functions ``sqrt, exp, log, log2, log10, sin, cos, tan, asin, acos, atan, sinh, cosh, tanh, asinh, acosh, atanh, acosh, atanh, erf, erfc, abs``.
simgi.cabal view
@@ -1,5 +1,5 @@ Name: simgi-Version: 0.1.1+Version: 0.2 License: GPL license-file: COPYING copyright: (C) 2009 Markus Dittrich@@ -9,16 +9,16 @@ molecular systems using Gillespie's method. Author: <haskelladdict@users.sourceforge.net> Maintainer: <haskelladdict@users.sourceforge.net>-stability: alpha+stability: beta build-type: Simple Homepage: http://simgi.sourceforge.net/ cabal-version: >= 1.6 extra-source-files: README Executable simgi- Build-Depends: base, containers >= 0.1.0.0,+ Build-Depends: base >= 2 && <= 4, containers >= 0.1.0.0, parsec == 2.1.*, mtl >= 1.1.0.2, haskell98,- random >= 1.0.0.1+ random >= 1.0.0.1, mersenne-random-pure64 >= 0.2 ghc-options: -O2 Main-Is: simgi.hs hs-source-dirs: src
src/CommandLine.hs view
@@ -20,66 +20,101 @@ -- | main gsim driver module CommandLine ( process_commandline - , Options(..)+ , SimgiOpt(..) ) where + -- imports+import Data.Word import Prelude import System import System.Console.GetOpt -- local imports+import GenericModel import Messages + -- | main driver for command line processing-process_commandline :: [String] -> IO ((CmdlRequest,String), [String])-process_commandline args = +process_commandline :: ModelState -> [String] + -> IO (ModelState, [String])+process_commandline state args = let (actions, nonOpts, _) = getOpt RequireOrder options args in- foldl (>>=) ( return defaultOptions ) actions >>= \opts ->+ foldl (>>=) ( return [] ) actions >>= \opts -> let - Options { cmdlRequest = request - , cmdlString = pattern- } = opts+ newState = process_options opts state in- return ((request,pattern),nonOpts) + return (newState,nonOpts) --- | possible options for commandline-data CmdlRequest = None | Help+-- | process all user provided commandline options and adjust+-- the model state accordingly+process_options :: [SimgiOpt] -> ModelState -> ModelState+process_options [] state = state+process_options ( (SimgiOpt { cmdlRequest = req+ , cmdlString = val }):xs) state = + case req of + Seed -> process_options xs ( state { seed = parse_seed val } )+ _ -> process_options xs state -- ignore unknown requests --- | data structure for keeping track of --- selected command line options-data Options = Options {+ + where+ parse_seed aSeed = floor (read aSeed :: Double) :: Word64 ++++-- | data structure describing all commandline options we know of+data CmdlRequest = Seed ++++-- | data structure for keeping track of a specific user provided+-- command line switch +data SimgiOpt = SimgiOpt { cmdlRequest :: CmdlRequest, cmdlString :: String } --- | default selections-defaultOptions :: Options-defaultOptions = Options {- cmdlRequest = None,- cmdlString = ""-} - -- | available command line flags-options :: [OptDescr (Options -> IO Options)]+options :: [OptDescr ([SimgiOpt] -> IO [SimgiOpt])] options = [ Option ['v'] ["version-info"] (NoArg version_info) - "show version information"]+ "show version information",+ Option ['h'] ["help"] (NoArg help_msg) "show help message",+ Option ['s'] ["seed"] (ReqArg seed_value "SEED") "seed value"+ ] + -- | extractor function for version info-version_info :: Options -> IO Options+version_info :: [SimgiOpt] -> IO [SimgiOpt] version_info _ = do show_version exitWith ExitSuccess++++-- | extractor function for help message+help_msg :: [SimgiOpt] -> IO [SimgiOpt]+help_msg _ =+ do+ usage+ exitWith ExitSuccess++++-- | extract the seed value +seed_value :: String -> [SimgiOpt] -> IO [SimgiOpt]+seed_value arg opt = + return ( (SimgiOpt { cmdlString = arg, cmdlRequest = Seed }):opt)+
src/Engine.hs view
@@ -19,25 +19,30 @@ --------------------------------------------------------------------} -- | the main compute Engine-module Engine ( create_initial_output+module Engine ( compute_trigger+ , create_initial_output , create_initial_state+ , execute_actions , gillespie_driver , module GenericModel ) where + -- imports import Control.Monad.State import qualified Data.Map as M import Prelude import Text.Printf+import System.Random(randomR)+import qualified System.Random.Mersenne.Pure64 as MT import System.IO + -- local imports import ExtraFunctions import GenericModel import RpnCalc --- import Debug.Trace -- | main simulation driver -- the simulator either stops when@@ -46,59 +51,82 @@ -- zero t_max is treated as being infinity gillespie_driver :: Handle -> Double -> Integer -> ModelState -> IO () gillespie_driver handle simTime dmpIter state = - let (output, outState) = runState run_gillespie $ state - (curTime, newState) = update_state dmpIter outState+ let + (output, outState) = runState run_gillespie $ state + (curTime, newState) = update_state dmpIter outState+ reversedOutput = reverse output in- (write_output handle . reverse $ output)- >> if curTime >= simTime- then return ()- else gillespie_driver handle simTime dmpIter newState+ -- write output to console and the output file+ (write_info $ head reversedOutput)+ >> (write_data handle reversedOutput) + -- next iteration if we're not at the end+ >> if curTime < simTime+ then gillespie_driver handle simTime dmpIter newState+ else return () + -- | updates the state for the next iteration update_state :: Integer -> ModelState -> (Double,ModelState) -update_state dataDumpIter state@(ModelState { currentTime = t - , maxIter = it- }) =- (t, state { maxIter = it + dataDumpIter, outputList = [] })+update_state dataDumpIter + state@(ModelState { currentTime = t + , outputBufferSize = it+ }) + = (t, state { outputBufferSize = it + dataDumpIter, outputCache = [] }) + -- | actual compute loop run_gillespie :: GillespieState [Output] run_gillespie = get - >>= \inState@(ModelState { molCount = in_mols- , reactions = in_reacts- , randNums = (r1:r2:randRest)- , currentTime = t- , currentIter = it- , maxTime = t_max- , maxIter = it_max- , outputFreq = freq- , outputList = output+ >>= \inState@(ModelState { molCount = in_mols+ , reactions = in_reacts+ , randGen = rGen+ , events = molEvents+ , currentTime = t+ , currentIter = it+ , maxTime = t_max+ , outputBufferSize = it_max+ , outputFreq = freq+ , outputRequest = outputVars+ , outputCache = output+ , variables = theVars }) -> + -- compute and update the next state- let out_rates = compute_rates in_reacts in_mols t []- a_0 = sum out_rates- tau = (-1.0/a_0) * log(r1)- t_new = t+tau- mu = get_mu (a_0*r2) out_rates- out_mols = adjust_molcount in_mols in_reacts mu- new_output = generate_output freq it t_new out_mols output- newState = inState { molCount = out_mols- , rates = out_rates- , randNums = randRest- , currentTime = t_new- , currentIter = it+1- , outputList = new_output- }+ let + -- generate two random numbers+ (r1,rGen1) = randomR (0.0 :: Double, 1.0) rGen+ (r2,rGen2) = randomR (0.0 :: Double, 1.0) rGen1++ -- update state+ symbols = SymbolTable in_mols theVars+ out_rates = compute_rates symbols in_reacts t []+ a_0 = sum out_rates+ tau = (-1.0/a_0) * log(r1)+ t_new = t+tau+ mu = get_mu (a_0*r2) out_rates+ out_mols = update_molcount in_mols in_reacts mu+ newSymbols = (symbols { molSymbols = out_mols })+ evt_syms = handle_events molEvents newSymbols t_new+ new_output = generate_output freq it t_new newSymbols outputVars output+ newState = inState { molCount = (molSymbols evt_syms)+ , rates = out_rates+ , randGen = rGen2+ , currentTime = t_new+ , currentIter = it+1+ , outputCache = new_output+ , variables = (varSymbols evt_syms)+ } in -- this prevents simulation from getting stuck -- FIXME: We need to come up with mechanism to propagate- -- error message corresponding to cases such as this one to the user!+ -- error message corresponding to cases such as this one + -- to the user! if (is_equal tau 0.0) then let finalState = newState { currentTime = t_max } in put finalState >> return output@@ -108,103 +136,232 @@ else put newState >> run_gillespie ++-- | handle all user defined events and return the adjusted+-- number of molecules+-- WARNING: We should probably check the Event Stack before we use+-- it to compute stuff; at least make sure molecule exist+handle_events :: [Event] -> SymbolTable -> Double -> SymbolTable+handle_events [] symbols _ = symbols+handle_events (x:xs) symbols t = + let+ newSymbols = handle_single_event x symbols t+ in+ handle_events xs newSymbols t++++-- | handle a single user event+handle_single_event :: Event -> SymbolTable -> Double -> SymbolTable+handle_single_event evt symbols t =+ let+ triggers = evtTrigger evt+ actions = evtActions evt+ triggerVal = compute_trigger symbols t triggers+ in + if triggerVal+ then execute_actions actions symbols t + else symbols++++-- | compute the value of a trigger+compute_trigger :: SymbolTable -> Double + -> ([EventTriggerPrimitive],[EventTriggerCombinator]) -> Bool+compute_trigger _ _ ([],_) = False -- this is should never happen+compute_trigger symbols t ((x:xs),combs) = compute_trigger_h (eval_trigger x) xs combs++ where+ compute_trigger_h acc [] _ = acc+ compute_trigger_h acc _ [] = acc+ compute_trigger_h acc (y:ys) (c:cs) = ++ case c of+ AndCombinator -> compute_trigger_h (acc && (eval_trigger y)) ys cs+ OrCombinator -> compute_trigger_h (acc || (eval_trigger y)) ys cs++ + eval_trigger e = (trigRelation e) (leftTrigger e) (rightTrigger e)+ leftTrigger = rpn_compute symbols t . trigLeftExpr + rightTrigger = rpn_compute symbols t . trigRightExpr++++-- | handle all actions associated with a user event+execute_actions :: [EventAction] -> SymbolTable -> Double + -> SymbolTable+execute_actions [] symbols _ = symbols+execute_actions (x:xs) symbols t =+ let+ newSymbol = execute_single_action x symbols t+ in+ execute_actions xs newSymbol t++++-- | handle a single event triggered action+execute_single_action :: EventAction -> SymbolTable -> Double+ -> SymbolTable+execute_single_action eventAction symbols t =+ let+ aName = evtName eventAction+ action = evtAct eventAction+ in+ case action of+ Constant c -> adjust_count aName c ++ Function rpn -> let+ newCount = rpn_compute symbols t rpn+ in+ adjust_count aName newCount + + where+ -- NOTE: presently, converting double -> int is done+ -- via floor. Is this a good policy (once documented+ -- properly)?+ to_int :: Double -> Int+ to_int = floor ++ -- adjust either a molecule count or the value of a+ -- variable+ adjust_count key val = case M.member key (molSymbols symbols) of++ True -> symbols { molSymbols = + M.insert key (to_int val) (molSymbols symbols) }+ False -> symbols { varSymbols = + M.insert key (Constant val) (varSymbols symbols) }+++ -- | generate a new Output data structure based on the current -- molecule counts-generate_output :: Integer -> Integer -> Double -> MoleculeMap - -> [Output] -> [Output]-generate_output afreq it t amol outlist +generate_output :: Integer -> Integer -> Double -> SymbolTable+ -> [String] -> [Output] -> [Output]+generate_output afreq it t symTable outVars outlist | mod it afreq /= 0 = outlist- | otherwise = new_out:outlist+ | otherwise = new_out:outlist where- new_out = Output { iteration = it- , time = t- , mols = amol+ currentOutputList = grab_output_data outVars t symTable++ new_out = Output { iteration = it+ , time = t+ , outputData = currentOutputList } +-- | given a list of variable or molecule names, goes through the+-- symbol table, grabs the current values associated with the variables,+-- and returns them as a list+grab_output_data :: [String] -> Double -> SymbolTable -> [Double]+grab_output_data vars aTime symbols = + foldr (\x acc -> (get_val_from_symbolTable x aTime symbols):acc) [] vars+ ++ -- | depending on which reaction happened adjust the number of -- molecules in the system-adjust_molcount :: MoleculeMap -> [Reaction] -> Int -> MoleculeMap-adjust_molcount theMap rs mID =+update_molcount :: MoleculeMap -> [Reaction] -> Int -> MoleculeMap+update_molcount theMap rs mID = - let (Reaction { react = react_in }) = rs !! mID+ let + (Reaction { reaction = react_in }) = rs !! mID in adjustMap react_in theMap where adjustMap :: [(String,Int)] -> MoleculeMap -> MoleculeMap adjustMap [] m = m- adjustMap ((k,a):rands) m = let val = (M.!) m k+ adjustMap ((k,a):changes) m = let + val = (M.!) m k m_new = M.insert k (a+val) m- in- adjustMap rands m_new+ in+ adjustMap changes m_new + -- | pick the \mu value for the randomly selected next reaction -- reaction to happen get_mu :: Double -> [Double] -> Int get_mu val = length . takeWhile ( <val ) . scanl1 (+) - --- | compute the current value for the reaction probabilities based on --- the number of molecules and reaction rates-compute_rates :: [Reaction] -> MoleculeMap -> Double +++-- | compute the current value for the reaction probabilities based +-- on the number of molecules and reaction rates+compute_rates :: SymbolTable -> [Reaction] -> Double -> RateList -> RateList-compute_rates [] _ _ rts = reverse rts-compute_rates ((Reaction {rate = c_in, aList = a_in }):rs) - theMap theTime rts = +compute_rates _ [] _ rts = reverse rts+compute_rates symbols ((Reaction {rate = c_in, actors = a_in }):rs) + theTime rts = case c_in of- (Constant aRate) -> compute_rates rs theMap theTime+ (Constant aRate) -> compute_rates symbols rs theTime ((a_new aRate): rts)- (Function rateFunc) -> compute_rates rs theMap theTime- ((a_new . (rpn_compute theMap theTime) $ rateFunc):rts)+ (Function rateFunc) -> compute_rates symbols rs theTime+ ((a_new . (rpn_compute symbols theTime) $ rateFunc):rts) where mult = product $ map (\(a,f) -> f . fromIntegral $ - (M.!) theMap a) a_in + (M.!) (molSymbols symbols) a) a_in a_new = (*) mult + + -- | initialize the output data structure create_initial_output :: ModelState -> Output-create_initial_output (ModelState { molCount = initialMols }) = +create_initial_output (ModelState { molCount = initialMols + , variables = initialVars+ , outputRequest = outVars+ }) = Output { iteration = 1 , time = 0.0- , mols = initialMols+ , outputData = initialOutput } + where+ symbols = SymbolTable initialMols initialVars+ initialOutput = grab_output_data outVars 0.0 symbols+++ -- | set up the initial state-create_initial_state:: ModelState -> [Double] -> Output -> ModelState-create_initial_state state rand output = +create_initial_state:: ModelState -> Output -> ModelState+create_initial_state state@(ModelState { seed = theSeed}) out = state { rates = defaultRateList- , randNums = rand+ , randGen = MT.pureMT theSeed , currentTime = 0.0 , currentIter = 1- , outputList = [output]+ , outputCache = [out] } --- | basic routine writing the simulation output to stdout-write_output :: Handle -> [Output] -> IO ()-write_output _ [] = return ()-write_output handle ((Output {iteration = it, time = t, mols = m}):xs) = - let header = (printf "%-10d %18.15g" it t) :: String- counts = create_count_string m++-- | routine for writing basic accounting info to stdout+write_info :: Output -> IO ()+write_info (Output {iteration = it, time = t}) = + putStrLn $ printf "iteration: %-10d --> time: %6.5g s" it t ++++-- | basic routine writing the simulation output to the +-- file handle corresponding to the output file+write_data :: Handle -> [Output] -> IO ()+write_data _ [] = return ()+write_data handle ((Output {iteration = it, time = t, outputData = out}):xs) = ++ let + header = (printf "%-10d %18.15g " it t) :: String+ counts = create_count_string out in- -- write molecule data to output file hPutStrLn handle (header ++ counts)-- -- write current iteration to stdout- >> (putStrLn $ "iteration | time ---> " ++ header)- >> write_output handle xs+ >> write_data handle xs where- create_count_string :: MoleculeMap -> String- create_count_string = foldr (\x a -> (printf "%10d " x) ++ a) - "" . M.elems + create_count_string :: [Double] -> String+ create_count_string = foldr (\x a -> (printf "%18.15f " x) ++ a) ""
src/ExtraFunctions.hs view
@@ -26,9 +26,9 @@ , fact , is_equal , is_equal_with+ , maybe_to_int+ , maybe_to_positive_int , real_exp - , to_int- , to_positive_int ) where @@ -43,20 +43,25 @@ avogadroNum = 6.0221415e23 + -- | use glibc DBL_EPSILON dbl_epsilon :: Double dbl_epsilon = 2.2204460492503131e-16 ++ -- | comparison function for doubles via dbl_epsion is_equal :: Double -> Double -> Bool is_equal x y = abs(x-y) <= abs(x) * dbl_epsilon + -- | comparison function for doubles via threshold is_equal_with :: Double -> Double -> Double -> Bool is_equal_with x y th = abs(x-y) <= abs(x) * th + -- | function checking if a Double can be interpreted as a non -- negative Integer. We need this since all parsing of numbers -- is done with Doubles but some functions only work for @@ -65,22 +70,24 @@ -- Integer via floor and the compare if the numbers are identical. -- If yes, the number seems to be an Integer and we return it, -- otherwise Nothing-to_positive_int :: Double -> Maybe Integer-to_positive_int x = +maybe_to_positive_int :: Double -> Maybe Integer+maybe_to_positive_int x = case (is_equal (fromInteger . floor $ x) x) && (x > 0.0) of True -> Just $ floor x False -> Nothing + -- | function checking if a Double can be interpreted as an -- Integer. See is_positive_int for more detail-to_int :: Double -> Maybe Integer-to_int x = +maybe_to_int :: Double -> Maybe Integer+maybe_to_int x = case is_equal (fromInteger . floor $ x) x of True -> Just $ floor x False -> Nothing + -- | helper function for defining real powers -- NOTE: We use glibc's pow function since it is more -- precise than implementing it ourselves via, e.g.,@@ -92,10 +99,12 @@ real_exp a x = realToFrac $ c_pow (realToFrac a) (realToFrac x) + -- | factorial function fact :: Integer -> Integer fact 0 = 1 fact n = n * fact (n-1)+ -- | error function erf(x)
src/GenericModel.hs view
@@ -20,36 +20,85 @@ -- | data structures needed for defining a stochastic model module GenericModel ( defaultRateList+ , Event(..)+ , EventAction(..)+ , EventTriggerCombinator(..)+ , EventTriggerPrimitive(..) , GillespieState , initialModelState+ , MathExpr(..) , ModelState(..) , MoleculeMap , Output(..)- , Rate(..)+ , Rate , RateList , Reaction(..)+ , SymbolTable(..)+ , VariableMap+ , VariableValue ) where -- imports import Control.Monad.State import qualified Data.Map as M+import Data.List((\\))+import Data.Word import Prelude+import qualified System.Random.Mersenne.Pure64 as MT + -- local imports import RpnData +--import Debug.Trace ++ -- | A MoleculeMap keeps track of the current number of molecules type MoleculeMap = M.Map String Int --- | data type for reaction rates which can either be a Double--- or an RpnStack describing a function to compute the rate--- at run time-data Rate = Constant Double | Function RpnStack+-- | A VariableMap holds all definied variables and their +-- current value.+-- NOTE: variables may change their each iteration since+-- they may be time dependent.+type VariableMap = M.Map String MathExpr ++-- | SymbolTable holds all names we know about such as molecule+-- names, variable names ...+data SymbolTable = SymbolTable { molSymbols :: MoleculeMap+ , varSymbols :: VariableMap+ }+++-- | generic data type for a mathematical expression. This could+-- either be a constant or an expression inside an RpnStack+data MathExpr = Constant Double | Function RpnStack+++-- | make MathExpr an instance of Eq+-- We allow only comparison of constants with each other+-- and RPN stacks with each other+instance Eq MathExpr where+ + (Constant x) == (Constant y) = x == y+ (Function f1) == (Function f2) = f1 == f2+ _ == _ = False+++-- | data type for variable values which are of type MathExpr+-- i.e. they can either be a number or a function involving, e.g.,+-- TIME or molecule counts+type VariableValue = Double+++-- | data type for reaction rates which are of type MathExpr+type Rate = MathExpr++ -- | List of reactions and corresponding rates type RateList = [Double] @@ -57,63 +106,173 @@ defaultRateList = [] --- | for each elementary reaction i we need to provide --- 1) the reaction rate c_i or rate function --- 2) the reaction order (first, second, ...)--- 2) aList describing which molecular species are participating--- in a reaction (needed for computing h_mu in Gillespie's --- notation) and a function mapping a molecule count to the--- proper h_mu value (needed e.g. for cases where we have--- 2X terms where h_my would be 0.5*X*(X-1).--- 3) a list of tuple (i,j) describing that the count of molecule--- i changes by j should this reaction take place+-- | an actor is a description of a molecular species participating+-- in a reaction (needed for computing h_mu in Gillespie's +-- notation) and a function mapping a molecule count to the+-- proper h_mu value (needed e.g. for cases where we have+-- 2X terms where h_mu would be 0.5*X*(X-1).+type Actor = (String, Double -> Double)++++-- | for each elementary reaction i we need to keep track of+-- +-- rate : the reaction rate c_i or rate function +-- actors: a list of Actors+-- react : a list of tuples (i,j) describing that the reaction+-- changes the count of molecule i by j +-- data Reaction = Reaction { rate :: Rate- , aList :: [(String,Double -> Double)]- , react :: [(String,Int)]+ , actors :: [Actor]+ , reaction :: [(String,Int)] } +-- | make Reaction an instance of Eq so we can compare them+-- (used in our unit tests)+instance Eq Reaction where++ x == y = compare_reactions x y++ where+ -- | compare two reactions and return True if they+ -- are equal and false otherwise+ compare_reactions :: Reaction -> Reaction -> Bool+ compare_reactions + (Reaction { rate = rate1+ , actors = actors1+ , reaction = reaction1 })+ (Reaction { rate = rate2+ , actors = actors2 + , reaction = reaction2 }) =++ let+ rateComp = rate1 == rate2+ actorComp = compare_actors actors1 actors2++ -- since reaction is assembled from Data.Map in the+ -- input parser we can't use == here+ reactComp = ( reaction1 \\ reaction2 ) == []+ in+ rateComp && actorComp && reactComp+++ -- | compare two actors+ compare_actors :: [Actor] -> [Actor] -> Bool+ compare_actors [] [] = True+ compare_actors xs ys = and $ zipWith compare_actor_elem xs ys+++ -- | compare an two actor elements+ compare_actor_elem :: Actor -> Actor -> Bool+ compare_actor_elem act1 act2 =+ + let+ testNum = 133.0 :: Double+ nameComp = fst act1 == fst act2+ funcComp = (snd act1) testNum == (snd act2) testNum+ in+ nameComp && funcComp+++++-- | data type describing an action triggered during an event+-- It consists of a String tracking the molecule affected+-- as well as a mathematic expression describing the new molecule+-- count for this molecule+data EventAction = EventAction { evtName :: String+ , evtAct :: MathExpr+ }++++-- | data type describing an expression that triggers a +-- user event+data EventTriggerPrimitive = EventTriggerPrimitive + { trigLeftExpr :: RpnStack+ , trigRelation :: Double -> Double -> Bool+ , trigRightExpr :: RpnStack+ }+++-- | combinators that can be used to combine EventTriggerPrimitives+data EventTriggerCombinator = AndCombinator | OrCombinator+++-- | data type keeping track of possible events occuring during+-- the simulation. Each event consist of a+--+-- <trigger>: list of expressions each evaluating to a bool.+-- event is triggered if all expression evaluate to true+-- FIXME: In the future we should support more complex +-- boolen operations involving &&, ||, etc.+-- +-- <action>: a list of semicolon separated expressions of the form +--+-- mol/var = <numerical expression>+--+-- changing the value of mol/var by <numerical expression>+--+data Event = Event { evtTrigger :: ([EventTriggerPrimitive], [EventTriggerCombinator])+ , evtActions :: [EventAction]+ }+++ -- | Our model state-data ModelState = ModelState { molCount :: MoleculeMap- , rates :: RateList- , reactions :: [Reaction]- , randNums :: [Double]- , systemVol :: Double- , currentTime :: Double- , currentIter :: Integer- , maxTime :: Double- , maxIter :: Integer- , outputFreq :: Integer- , outputList :: [Output]- , outfileName :: String+data ModelState = ModelState { molCount :: MoleculeMap+ , rates :: RateList+ , reactions :: [Reaction]+ , randNums :: [Double]+ , seed :: Word64+ , randGen :: MT.PureMT+ , events :: [Event]+ , systemVol :: Double+ , currentTime :: Double+ , currentIter :: Integer+ , maxTime :: Double+ , outputBufferSize :: Integer+ , outputFreq :: Integer+ , outputRequest :: [String]+ , outputCache :: [Output]+ , outfileName :: String+ , variables :: VariableMap } type GillespieState a = State ModelState a + -- | data structure for keeping track of our output-data Output = Output { iteration :: Integer- , time :: Double- , mols :: MoleculeMap +data Output = Output { iteration :: Integer+ , time :: Double+ , outputData :: [Double] } deriving(Show) + -- | initial model state to be partially filled by the -- parser from the input deck initialModelState :: ModelState-initialModelState = ModelState { molCount = M.empty- , rates = []- , reactions = []- , randNums = []- , systemVol = 1.0- , currentTime = 0.0- , currentIter = 0- , maxTime = 0.0- , maxIter = 10000- , outputFreq = 1000- , outputList = []- , outfileName = ""+initialModelState = ModelState { molCount = M.empty+ , rates = []+ , reactions = []+ , randNums = []+ , events = []+ , seed = 1+ , randGen = MT.pureMT 1+ , systemVol = 1.0+ , currentTime = 0.0+ , currentIter = 0+ , maxTime = 0.0+ , outputBufferSize = 10000+ , outputFreq = 1000+ , outputRequest = []+ , outputCache = []+ , outfileName = ""+ , variables = M.empty }
src/InputCheck.hs view
@@ -37,45 +37,82 @@ -- should this ever become more extensive we should probably consider -- using Control.Monad.Error check_input :: ModelState -> Either String Bool -check_input (ModelState { molCount = theMols- , reactions = theReactions- , maxIter = iterCount- , outputFreq = outFreq- , outfileName = fileName+check_input (ModelState { molCount = theMols+ , reactions = theReactions+ , outputBufferSize = iterCount+ , outputFreq = outFreq+ , outfileName = fileName+ , variables = theVars+ , events = theEvents }) - = check_molecules (defined_mols theMols) - (react_mols theReactions)+ = check_molecules (M.keys theMols) (react_mols theReactions) >> check_positive_outfreq outFreq >> check_positive_itercount iterCount >> check_filename fileName- >> check_reaction_rate_functions (defined_mols theMols) - (rate_mols theReactions)+ >> check_variable_names defined_names+ (extract_variable_names (M.elems theVars)) "variables"+ >> check_variable_names defined_names + (extract_variable_names_from_rates theReactions) "reactions"+ >> check_variable_names defined_names+ (extract_variable_names_from_events theEvents) "events" where -- | extract all reaction participants- react_mols = L.nub . L.concat . map (map (fst) . react) + react_mols = L.nub . L.concat . map (map (fst) . reaction) - -- | extract all definied molecules- defined_mols = M.keys + -- | extract all definied names (molecules, variables, ...)+ defined_names = (M.keys theMols) ++ (M.keys theVars) + - -- | extract all molecules appearing in reaction rate- -- functions- rate_mols theRates = - let- stacks = foldr extract_rate_func [] . map rate $ theRates- in- L.nub . concat . map (foldr extract_rate_vars []) $ stacks+-- | extract all molecules/variables appearing in reaction rate+-- functions+extract_variable_names_from_rates :: [Reaction] -> [String]+extract_variable_names_from_rates = extract_variable_names . map rate - where- extract_rate_func (Function a) acc = a:acc- extract_rate_func _ acc = acc - extract_rate_vars (Variable a) acc = a:acc- extract_rate_vars _ acc = acc +-- | extract all variable/molecule names appearing in a list of+-- rpn stacks corresponding to some rate of variable definition+-- expression +extract_variable_names :: [MathExpr] -> [String]+extract_variable_names inputList =+ let+ stacks = foldr extract_rate_func [] inputList+ in+ L.nub . concat . map (foldr extract_rate_vars []) $ stacks + where+ extract_rate_func (Function a) acc = (toList a):acc+ extract_rate_func _ acc = acc++ extract_rate_vars (Variable a) acc = a:acc+ extract_rate_vars _ acc = acc++++-- | extract all variable/molecule names from expressions inside+-- events, i.e. insider triggers and actions+extract_variable_names_from_events :: [Event] -> [String]+extract_variable_names_from_events theEvents = ++ L.nub $ allActionNames theEvents +++ (extract_variable_names $ triggerExps theEvents ++ actionExps theEvents)++ where+ allTriggers = concat . foldr ((:) . fst . evtTrigger) []+ triggerExps = foldr (\x acc -> + ((Function $ trigLeftExpr x):(Function $ trigRightExpr x):acc)) + [] . allTriggers+ + allActions = concat . foldr ((:) . evtActions) []+ actionExps = foldr ((:) . evtAct) [] . allActions++ allActionNames = foldr ((:) . evtName) [] . allActions+++ -- | make sure the user specified an output file name check_filename :: String -> Either String Bool check_filename name @@ -83,6 +120,7 @@ | otherwise = Right True + -- | make sure all molecules in reactions are defined check_molecules :: [String] -> [String] -> Either String Bool check_molecules defMols reactMols = @@ -114,15 +152,15 @@ -- | make sure the user defined reaction rate function reference -- only existing molecule names-check_reaction_rate_functions :: [String] -> [String] - -> Either String Bool-check_reaction_rate_functions defMols rateMols =+check_variable_names :: [String] -> [String] -> String -> Either String Bool+check_variable_names defMols rateMols checkType = let - no_mol = rateMols L.\\ defMols+ noMol = rateMols L.\\ defMols in- case null no_mol of+ case null noMol of True -> Right True False -> Left $- "Error: The following molecules defined in reaction "- ++ "rates do not exist: " - ++ (L.concat $ L.intersperse "," no_mol)+ "Error: The following molecules or variables defined in the "+ ++ checkType+ ++ " block do not exist:\n --> " + ++ (L.concat $ L.intersperse ", " noMol)
src/InputParser.hs view
@@ -19,7 +19,10 @@ --------------------------------------------------------------------} -- | input file parser -module InputParser ( input_parser ) where+module InputParser ( input_parser+ , parse_events + , parse_reaction+ ) where -- imports import Control.Monad@@ -40,34 +43,209 @@ -- | main parser entry point input_parser :: CharParser ModelState ModelState input_parser = whiteSpace - *> many ( choice [ try parse_parameter_def+{- *> many ( choice [ try parse_parameter_def , try parse_molecule_def , try parse_reaction_def- ])+ , try parse_event_def+ ])-}+ *> parse_parameter_def+ *> optional (try parse_variable_def)+ *> parse_molecule_def+ *> parse_reaction_def+ *> optional (try parse_event_def)+ *> optional parse_output_def *> eof >> getState <?> "main parser" ++-- | parser for variable definitions+parse_variable_def :: CharParser ModelState ()+parse_variable_def = join ( updateState <$> insert_variables <$>+ parse_def_block "variables" (many parse_variable) ) + <?> "variable definition block" ++ where+ insert_variables :: [(String, MathExpr)] -> ModelState -> ModelState+ insert_variables theVars state = state { variables = M.fromList theVars }+++-- | parser for a single variable definition+parse_variable :: CharParser ModelState (String, MathExpr)+parse_variable = tuple_it <$> ((try parse_variable_name) )+ <*> (symbol "=" *> parse_variable_definition)+ <?> "variable definition"+ + where+ tuple_it one two = (one, two)+++-- | parser for variable name+parse_variable_name :: CharParser ModelState String+parse_variable_name = identifier + <?> "variable name"+++-- | parse the definition for a variable+parse_variable_definition :: CharParser ModelState MathExpr+parse_variable_definition = (try parse_constant_expression)+ <|> (braces parse_function_expression)+ <?> "variable value"+++-- | parser for output definitions+parse_output_def :: CharParser ModelState ()+parse_output_def = join ( updateState <$> insert_output_request <$>+ parse_def_block "output" (parse_output_list) ) + <?> "event definitions" ++ where+ insert_output_request :: [String] -> ModelState -> ModelState+ insert_output_request outDataList state = state { outputRequest = outDataList }+++-- | parse the list with variables or molecules to be punched to the +-- output file+parse_output_list :: CharParser ModelState [String]+parse_output_list = brackets (commaSep parse_variable_name)+++-- | parser for event definitions+parse_event_def :: CharParser ModelState ()+parse_event_def = join ( updateState <$> insert_events <$>+ parse_def_block "events" (many parse_events) ) + <?> "event definitions" ++ where+ insert_events :: [Event] -> ModelState -> ModelState+ insert_events newEvents state = state { events = newEvents }++++-- | parser for individual events+parse_events :: CharParser ModelState Event+parse_events = Event <$> (parse_trigger) <*> (reservedOp "=>" *> parse_actions)+ <?> "reaction event"++++-- | parser for an event trigger+parse_trigger :: CharParser ModelState + ([EventTriggerPrimitive], [EventTriggerCombinator])+parse_trigger = braces parse_trigger_expressions+ <?> "event trigger block"++++-- | parse a list of trigger expressions+parse_trigger_expressions :: CharParser ModelState + ([EventTriggerPrimitive], [EventTriggerCombinator])+parse_trigger_expressions = combine_it <$> parse_single_trigger_expression + <*> (many parse_boolean_trigger_expression)+ <?> "event trigger"+ + where+ combine_it e = foldr (\(x,y) (u,v) -> (x:u,y:v) ) ([e],[])++++-- | parse a single trigger expression+parse_single_trigger_expression :: CharParser ModelState EventTriggerPrimitive+parse_single_trigger_expression = + EventTriggerPrimitive <$> parse_infix_to_rpn <*> parse_relational+ <*> parse_infix_to_rpn+ <?> "single event trigger expression"++++-- | parse a single trigger expression prefixed with a && or ||+parse_boolean_trigger_expression :: CharParser ModelState + (EventTriggerPrimitive, EventTriggerCombinator)+parse_boolean_trigger_expression = + tuple_it <$> parse_boolean_combinator <*> parse_single_trigger_expression+ <?> "boolean trigger expression"++ where+ tuple_it a b = (b,a)++++-- | parse a boolean combinator (&& or ||)+parse_boolean_combinator :: CharParser ModelState EventTriggerCombinator+parse_boolean_combinator = try parse_AND <|> parse_OR+ <?> "boolean combinator"+ +++-- | parse an && combinator+parse_AND :: CharParser ModelState EventTriggerCombinator+parse_AND = symbol "&&" *> (pure AndCombinator)+ <?> "&&"+ +++-- | parse an || combinator+parse_OR :: CharParser ModelState EventTriggerCombinator+parse_OR = symbol "||" *> (pure OrCombinator)+ <?> "||"++++-- | parse a relational expression and return its associated+-- binary function+parse_relational :: CharParser ModelState (Double -> Double -> Bool)+parse_relational = try ( reservedOp ">=" >> pure (>=) )+ <|> try ( reservedOp "<=" >> pure (<=) )+ <|> try ( reservedOp "==" >> pure (==) )+ <|> ( reservedOp ">" >> pure (>) )+ <|> ( reservedOp "<" >> pure (<) )+ <?> "relational expression"++++-- | parser for an event action+parse_actions :: CharParser ModelState [EventAction]+parse_actions = braces parse_action_expressions+ <?> "event action block"++++-- | parser for a list of action expressions+parse_action_expressions :: CharParser ModelState [EventAction]+parse_action_expressions = + parse_single_action_expression `sepEndBy` semi + <?> "event action expression"+++-- | parser for a single event action expression+parse_single_action_expression :: CharParser ModelState EventAction+parse_single_action_expression = EventAction <$> + (molname) <*> (reservedOp "=" *> parse_function_expression)+ <?> "single event action expression"+++ -- | parser for simulation parameters parse_parameter_def :: CharParser ModelState ()-parse_parameter_def = between (reserved "def" *> reserved "parameters")- (reserved "end")- (parse_parameters `sepBy` whiteSpace) +parse_parameter_def = parse_def_block "parameters" (many parse_parameters) *> pure () <?> "parameter definitions" + -- | parse the individual parameters parse_parameters :: CharParser ModelState () parse_parameters = parse_time <|> parse_outputFile- <|> parse_outputIter+ <|> parse_outputBuffer <|> parse_outputFreq <|> parse_systemVol- <?> "time, outputIter, systemVol, outputFreq, outputFile"+ <?> "time, outputBuffer, systemVol, outputFreq,\+ \outputFile" + -- | parse the simulation time specs parse_time :: CharParser ModelState () parse_time = join (updateState <$> insert_time @@ -78,6 +256,7 @@ insert_time t state = state { maxTime = t } + -- | parse the value of the simulated system volume parse_systemVol :: CharParser ModelState () parse_systemVol = join (updateState <$> insert_volume@@ -94,6 +273,7 @@ insert_volume vol state = state { systemVol = vol } + -- | parse the name of the output file -- accepts paths but will NOT create any of the parents parse_outputFile :: CharParser ModelState ()@@ -102,25 +282,28 @@ *> parse_filename )) where- insert_filename name state = state { outfileName = name }+ insert_filename aName state = state { outfileName = aName } + -- | parse a filename parse_filename :: CharParser ModelState String parse_filename = stringLiteral + -- | parse the output iteration specification if present-parse_outputIter :: CharParser ModelState ()-parse_outputIter = join (updateState <$> insert_outputIter- <$> (reserved "outputIter" *> reservedOp "="+parse_outputBuffer :: CharParser ModelState ()+parse_outputBuffer = join (updateState <$> insert_outputBuffer+ <$> (reserved "outputBuffer" *> reservedOp "=" *> integer )) where- insert_outputIter i state = state { maxIter = i }+ insert_outputBuffer i state = state { outputBufferSize = i } + -- | parse the output iteration specification if present parse_outputFreq :: CharParser ModelState () parse_outputFreq = join (updateState <$> insert_outputFreq@@ -131,12 +314,11 @@ insert_outputFreq i state = state { outputFreq = i } + -- | parser for molecule definitions parse_molecule_def :: CharParser ModelState () parse_molecule_def = join ( updateState <$> insert_molecules <$> - between (reserved "def" *> reserved "molecules")- (reserved "end")- (parse_molecules `sepBy` whiteSpace) )+ parse_def_block "molecules" (many parse_molecules)) <?> "molecule definitions" where@@ -145,9 +327,11 @@ state { molCount = M.fromList theMols } + -- | parse a molecule name and the number of molecules of this type parse_molecules :: CharParser ModelState (String,Int)-parse_molecules = make_molecule <$> (try molname) <*> integer+parse_molecules = make_molecule <$> (try molname) <*> (symbol "=" *> integer)+ <?> "molecule expression" where make_molecule mol aCount = (mol,fromInteger aCount) @@ -156,24 +340,14 @@ -- A molecule name can consist of letters and numbers but has to -- start with a letter. The following keywords are reserved molname :: CharParser ModelState String-molname = not_end ((:) <$> letter <*> many (alphaNum <?> ""))- <?> "molecule name" ----- | short checker making sure we don't scan beyond the "end" statement--- of a block -not_end :: CharParser ModelState String -> CharParser ModelState String-not_end p = p >>= \name -> case name /= "end" of- True -> pure name- False -> pzero+molname = identifier+ <?> "molecule name" -- | parser for reaction definitions parse_reaction_def :: CharParser ModelState () parse_reaction_def = join ( updateState <$> insert_reactions <$>- between (reserved "def" *> reserved "reactions")- (reserved "end")- (parse_reaction `sepBy` whiteSpace) )+ parse_def_block "reactions" (many parse_reaction) ) <?> "reaction definitions" where@@ -181,6 +355,7 @@ insert_reactions reacts state = state { reactions = reacts } + -- | parser for a single reaction specification of the type -- aA + bB + cC + .... -> n1P1 + n2P2 + ...... : rate : -- NOTE: In order to convert the reaction rates (if requested@@ -190,7 +365,7 @@ parse_reaction = setup_reaction <$> (parse_react_prod <* reservedOp "->") <*> parse_react_prod - <*> parse_rate+ <*> parse_rate_expression <*> (getState >>= \(ModelState {systemVol = vol}) -> pure vol) where@@ -199,13 +374,13 @@ let action = create_react r p hFactor = create_hFact r - theRate = if (vol < 0.0) -- no rate conversion for systemVol = nil- then cin+ theRate = if (vol < 0.0) -- no rate conversion for + then cin -- systemVol = nil else convert_rate cin (M.size r) vol in - Reaction { rate = theRate- , aList = hFactor- , react = action+ Reaction { rate = theRate+ , actors = hFactor+ , reaction = action } @@ -222,21 +397,25 @@ case order of 1 -> theFunc _ -> let mult = 1.0/(avogadroNum * volume^(order-1)) in- Function $ stack ++ [Number mult,BinFunc (*)]+ Function . RpnStack $ (toList stack) + ++ [Number mult,BinFunc (*)] -- | create the list holding the molecule number changes for - -- this reaction+ -- this reactioni. If the net change in molecule number is+ -- zero ( a nop) we remove the action. create_react r p = let reacts = M.map (*(-1)) r in- M.assocs $ M.unionWith (+) reacts p+ M.assocs . snd . M.partition (==0)+ $ M.unionWith (+) reacts p+ -- | create the list containing the h factors- -- WARNING/FIXME: Currently, things are ill defined if the number of- -- molecules for species A is below the stoichiometric reactant+ -- WARNING/FIXME: Currently, things are ill defined if the number + -- of molecules for species A is below the stoichiometric reactant -- coefficient; i.e. if #A = 2 then 3A -> ... does not make sense create_hFact :: (M.Map String Int) -> [(String, Double -> Double)] create_hFact = create_hFact_h [] . M.assocs @@ -256,34 +435,14 @@ generate_lambda n x = (x-n+1) * generate_lambda (n-1) x --- | parse rate parses a reaction rate. This can either be a simple--- constant of a full blown infix math expression.--- Reaction rates must be enclosed by colons ":"-parse_rate :: CharParser ModelState Rate-parse_rate = (try parse_constant_rate) <|> parse_rate_function- <?> "rate constant or rate function"----- | parser for a simple rate constant expression-parse_constant_rate :: CharParser ModelState Rate-parse_constant_rate = Constant <$> braces parse_number- <?> "rate constant" ----- | parser for a rate function-parse_rate_function :: CharParser ModelState Rate-parse_rate_function = Function <$> braces parse_infix_to_rpn- <?> "rate function" -- -- | parse list of reactants/products of reaction -- we expect to parse a stream that looks like -- n_1 R1 + n_2 R2 + n_3 R3 + ..... -- If n_i is missing we assume it is 1.0 parse_react_prod :: CharParser ModelState (M.Map String Int) parse_react_prod = (reserved "nil" *> pure (M.empty))- <|> (M.fromList <$> ((make_tuple <$> option 1 integer - <*> (try molname <* whiteSpace)) + <|> (M.fromList <$> + ((make_tuple <$> option 1 integer <*> try molname) `sepBy` reservedOp "+") ) <?> "reactant or product list" @@ -291,6 +450,7 @@ make_tuple x y = (y, fromInteger x) + -- | parse a number, can be used with 'many' and other parser -- combinators; integers are automatically promoted to double parse_number :: CharParser ModelState Double@@ -302,6 +462,7 @@ Right x -> x + -- | parse a positive number, can be used with 'many' and other -- parser combinators; integers are automatically promoted to double parse_positive_number :: CharParser ModelState Double@@ -315,3 +476,45 @@ else pzero <?> "unsigned integer or double" ++++-- | parser for a def block structure+parse_def_block :: String -> CharParser ModelState a + -> CharParser ModelState a+parse_def_block blockName parser = ++ between (reserved "def" *> reserved blockName )+ (reserved "end")+ (parser)+ <?> "parameter definitions"++++-- | parse a simple rate expression+-- FIXME: We can not re-use parse_expression below +-- since currently the order of function/constant parsing +-- has to be reversed otherwise rates are always parsed+-- as trivial Functions. +parse_rate_expression :: CharParser ModelState MathExpr+parse_rate_expression = (try (lineToken parse_constant_expression))+ <|> (lineToken (braces parse_function_expression)) + <?> "constant or function expression"++ where+ lineToken = between (symbol "|") (symbol "|") +++-- | parser for a simple rate constant expression+parse_constant_expression :: CharParser ModelState MathExpr+parse_constant_expression = Constant <$> parse_number+ <?> "rate constant" ++++-- | parser for a rate function+parse_function_expression :: CharParser ModelState MathExpr+parse_function_expression = Function <$> parse_infix_to_rpn+ <?> "rate function" ++
src/Messages.hs view
@@ -24,22 +24,57 @@ , usage ) where + -- imports import Prelude +-- local imports+import GenericModel++ -- | show version info show_version :: IO ()-show_version = putStrLn "This is simgi v0.1.1 (C) 2009 Markus Dittrich"+show_version = putStrLn "This is simgi v0.2 (C) 2009 Markus Dittrich" + -- | show a brief startup message-startup_message :: IO ()-startup_message = show_version - >> putStrLn "\nstarting simulation ..... here we go\n"+startup_message :: ModelState -> IO ()+startup_message state = show_version+ >> (putStrLn $ "\n-------- Simulation parameters ----------")+ >> (putStrLn $ "rng seed : " ++ (show simSeed))+ >> (putStrLn $ "max time : " ++ (show simTime) ++ " s")+ >> (putStrLn $ "system volume : " ++ finalSimVol)+ >> (putStrLn $ "data output frequency : " ++ (show simFreq))+ >> (putStrLn $ "log output frequency : " ++ (show simLogFreq))+ >> (putStrLn $ "output filename : " ++ simOutFile)+ >> (putStrLn $ "-----------------------------------------")+ >> putStrLn "\nstarting simulation ...\n" + where+ (ModelState { seed = simSeed + , maxTime = simTime+ , systemVol = simVol+ , outputBufferSize = simLogFreq+ , outputFreq = simFreq+ , outfileName = simOutFile + }) = state + -- for nil volumes we want to display nil not -1.0+ finalSimVol = if (simVol < 0) + then "nil"+ else (show simVol) ++ " m^3"++ -- | provide brief usage info usage :: IO ()-usage = putStrLn "Usage: simgi <input file>"+usage = putStrLn "Usage: simgi [options] <input file>\n\n\+ \Currently supported options are:\n\n\+ \\t -s --seed <seed value> \n\+ \\t specify value for the rng starting seed.\n\n\+ \\t -v --version-info \n\+ \\t print version info.\n\n\+ \\t -h --help \n\+ \\t print this help message."
src/RpnCalc.hs view
@@ -21,7 +21,9 @@ -- | RpnCalc defines the data structures and a calculator engine -- for computing mathematical expressions that have been parsed -- into reverse polish notations-module RpnCalc ( rpn_compute ) where+module RpnCalc ( rpn_compute+ , get_val_from_symbolTable + ) where -- imports @@ -32,14 +34,16 @@ import GenericModel import RpnData +-- import Debug.Trace -- | computes an expressions based on an rpn stack+-- molecule names are looked up in a MoleculeMap -- NOTE: This function expects the RPNstack to be sanitized -- with respect to the variables, i.e., all variables in--- the stack are assumed to exist in the MoleculeMap-rpn_compute :: MoleculeMap -> Double -> RpnStack -> Double-rpn_compute _ _ [(Number x)] = x-rpn_compute molMap theTime xs = num +-- the stack are assumed to exist in the VariableMap+rpn_compute :: SymbolTable -> Double -> RpnStack -> Double+rpn_compute _ _ (RpnStack [(Number x)]) = x+rpn_compute symbols theTime (RpnStack xs) = num where (Number num) = head . foldl evaluate [] $ xs@@ -56,9 +60,18 @@ evaluate ys (Time) = (Number theTime):ys -- extract molecule variable- evaluate ys (Variable x) = (Number $ replace_var x):ys- where- replace_var :: String -> Double- replace_var = fromIntegral . (M.!) molMap+ evaluate ys (Variable x) = (Number $ get_val_from_symbolTable x theTime symbols):ys evaluate ys item = item:ys+++-- | retrieve the value of a given symbol (either variable or molecule count) from+-- the symbol table +get_val_from_symbolTable :: String -> Double -> SymbolTable -> Double+get_val_from_symbolTable var aTime symbols =+ + case M.lookup var (molSymbols symbols) of+ Just value -> fromIntegral value+ Nothing -> case (M.!) (varSymbols symbols) var of+ Constant c -> c+ Function s -> rpn_compute symbols aTime s
src/RpnData.hs view
@@ -22,7 +22,7 @@ -- for computing mathematical expressions that have been parsed -- into reverse polish notations module RpnData ( RpnItem(..)- , RpnStack+ , RpnStack(..) ) where @@ -38,9 +38,58 @@ | UnaFunc (Double -> Double) | BinFunc (Double -> Double -> Double) -type RpnStack = [RpnItem] +-- | RpnStack describes a computation stored in a stack of+-- RpnItems+newtype RpnStack = RpnStack { toList :: [RpnItem] } ++-- | make RpnStack an instance of Eq. +-- We compare two RpnStacks by computing them and replacing+-- each occuring Variable by a default value. This should+-- be ok in all but some pathological cases.+instance Eq RpnStack where+ + x == y = (internal_rpn_compute x) == (internal_rpn_compute y)+++-- | computes an expressions based on an rpn stack+-- NOTE: This function is intended to be used for+-- the Eq instance of RpnItem only!+-- names are replaced by defaultVar. ++-- | default variable used to replace all encountered+-- variable. Its value is arbitrary+defaultVar :: Double+defaultVar = 47.0++internal_rpn_compute :: RpnStack -> Double+internal_rpn_compute (RpnStack [(Number x)]) = x+internal_rpn_compute (RpnStack xs) = num ++ where+ (Number num) = head . foldl evaluate [] $ xs++ -- evaluate unary function (sin, cos, ..)+ evaluate ((Number x):ys) (UnaFunc f) = + (Number $ f x):ys++ -- evaluate binary function (*,+,..)+ evaluate ((Number x):(Number y):ys) (BinFunc f) =+ (Number $ f y x):ys++ -- extrace current time+ evaluate ys (Time) = (Number defaultVar):ys++ -- extract molecule variable+ evaluate ys (Variable x) = (Number defaultVar):ys++ evaluate ys item = item:ys++++-- | a pretty lame show instance+-- use only for debugging purposes instance Show RpnItem where show (Time) = "TIME" show (Number x) = show x
src/RpnParser.hs view
@@ -37,12 +37,13 @@ -- | parses a mathematical infix expression and converts -- into a stack in rpn parse_infix_to_rpn :: CharParser ModelState RpnStack -parse_infix_to_rpn = add_term -- >>= \inp -> trace (show inp) (return inp)+parse_infix_to_rpn = RpnStack <$> add_term <?> "infix math expression" + -- | parser for expressions chained via "+" -add_term :: CharParser ModelState RpnStack+add_term :: CharParser ModelState [RpnItem] add_term = concat . insert_adds <$> sub_term `sepBy` (reservedOp "+") <?> "addition term"@@ -52,8 +53,9 @@ insert_adds (x:xs) = x:foldr (\y a -> y:[BinFunc (+)]:a) [] xs + -- | parser for expressions chained via "-"-sub_term :: CharParser ModelState RpnStack+sub_term :: CharParser ModelState [RpnItem] sub_term = concat . insert_subs <$> div_term `sepBy` (reservedOp "-") <?> "subtraction term"@@ -63,8 +65,9 @@ insert_subs (x:xs) = x:foldr (\y a -> y:[BinFunc (-)]:a) [] xs + -- | parser for expressions chained via "*" -div_term :: CharParser ModelState RpnStack+div_term :: CharParser ModelState [RpnItem] div_term = concat . insert_divs <$> mul_term `sepBy` (reservedOp "/") <?> "division term"@@ -74,8 +77,9 @@ insert_divs (x:xs) = x:foldr (\y a -> y:[BinFunc (/)]:a) [] xs + -- | parser for expressions chained via "/"-mul_term :: CharParser ModelState RpnStack+mul_term :: CharParser ModelState [RpnItem] mul_term = concat . insert_muls <$> exp_term `sepBy` (reservedOp "*") <?> "product term"@@ -85,9 +89,10 @@ insert_muls (x:xs) = x:foldr (\y a -> y:[BinFunc (*)]:a) [] xs + -- | parser for potentiation operations "^"-exp_term :: CharParser ModelState RpnStack-exp_term = concat . insert_exps <$> (whiteSpace *> factor) `sepBy` (reservedOp "^")+exp_term :: CharParser ModelState [RpnItem]+exp_term = concat . insert_exps <$> factor `sepBy` (reservedOp "^") <?> "exponent" where@@ -95,9 +100,10 @@ insert_exps (x:xs) = x:foldr (\y a -> y:[BinFunc real_exp]:a) [] xs + -- | parser for individual factors, i.e, numbers, -- variables or operations-factor :: CharParser ModelState RpnStack+factor :: CharParser ModelState [RpnItem] factor = try parse_single_number -- need try due to the unary "-" <|> try signed_parenthesis -- (otherwise we get stuck) <|> parse_functions@@ -105,9 +111,10 @@ <?> "token or variable" + -- | parse all operations of type (Double -> Double) -- we currently know about-parse_functions :: CharParser ModelState RpnStack+parse_functions :: CharParser ModelState [RpnItem] parse_functions = (msum $ extract_ops builtinFunctions) <?> "builtin unary function" where@@ -124,10 +131,11 @@ insert_func f xs = xs ++ [UnaFunc f] + -- | parse a potentially signed expression enclosed in parenthesis. -- In the case of parenthesised expressions we -- parse -(...) as (-1.0)*(...)-signed_parenthesis :: CharParser ModelState RpnStack+signed_parenthesis :: CharParser ModelState [RpnItem] signed_parenthesis = push_parens <$> parse_sign <*> parens add_term <?> "signed parenthesis" @@ -135,17 +143,19 @@ push_parens sign xs = xs ++ [Number sign,BinFunc (*)] + -- | parse a single number; integers are automatically promoted -- to double -- NOTE: Due to the notFollowedBy this parser can not be used -- with 'many' and other parser combinators.-parse_single_number :: CharParser ModelState RpnStack-parse_single_number = push <$> (parse_number <* notFollowedBy alphaNum)+parse_single_number :: CharParser ModelState [RpnItem]+parse_single_number = push <$> (parse_number) <?> "signed integer or double" where push x = [Number x] + -- | parse a number, can be used with 'many' and other parser -- combinators; integers are automatically promoted to double parse_number :: CharParser ModelState Double@@ -157,17 +167,17 @@ Right x -> sign * x + -- | parse the sign of a numerical expression parse_sign :: CharParser ModelState Double-parse_sign = option 1.0 ( whiteSpace *> char '-' *> pure (-1.0) )+parse_sign = option 1.0 ( char '-' *> pure (-1.0) ) <?> "sign" + -- | this is how valid variable names have to look like-parse_variable :: CharParser ModelState RpnStack-parse_variable = - push <$> parse_sign - <*> ((:) <$> letter <*> many (alphaNum <?> "") <* whiteSpace)+parse_variable :: CharParser ModelState [RpnItem]+parse_variable = push <$> parse_sign <*> identifier <?> "variable" where -- in case of a unary minus we also push the necessary
src/TokenParser.hs view
@@ -22,15 +22,20 @@ module TokenParser ( module Control.Applicative , module Text.ParserCombinators.Parsec , braces+ , brackets , builtinFunctions+ , colon , comma+ , commaSep , charLiteral , float+ , identifier , integer , parens , keywords , lexer , naturalOrFloat+ , operator , operators , reservedOp , reserved@@ -48,7 +53,6 @@ import Text.ParserCombinators.Parsec hiding (many,optional, (<|>)) import qualified Text.ParserCombinators.Parsec.Token as PT import Text.ParserCombinators.Parsec.Language (haskellDef- , opLetter , reservedOpNames , reservedNames ) @@ -64,6 +68,7 @@ (<*>) = ap + -- |Alternative instance for MonadPlus instance Alternative (GenParser s a) where empty = mzero@@ -98,36 +103,47 @@ ] + -- | all other keywords that are not regular functions keywords :: [String] keywords = [ "def", "molecules", "reactions", "time", "outputIter"- , "nil", "outputFreq", "outputFile", "systemVol"+ , "nil", "outputFreq", "outputFile", "systemVol", "seed"+ , "end", "variables", "time" ] operators :: [String]-operators = ["+","->","::","=","{","}"]+operators = ["+", "->", "::", "=", "{", "}", ">=", "==", "<="+ , "<", ">", "*", "/", "-", "&&"] + -- | function generating a token parser based on a -- lexical parser combined with a language record definition lexer :: PT.TokenParser st lexer = PT.makeTokenParser ( haskellDef { reservedOpNames = operators- , opLetter = oneOf "*+/^" , reservedNames = keywords ++ map fst builtinFunctions } ) + -- | token parser for parenthesis parens :: CharParser st a -> CharParser st a parens = PT.parens lexer ++ -- | token parser for parenthesis braces :: CharParser st a -> CharParser st a braces = PT.braces lexer +-- | token parser for brackets +brackets :: CharParser st a -> CharParser st a+brackets = PT.brackets lexer++ -- | token parser for Integer integer :: CharParser st Integer integer = PT.integer lexer@@ -158,27 +174,53 @@ reservedOp = PT.reservedOp lexer + -- | token parser for keywords reserved :: String -> CharParser st () reserved = PT.reserved lexer + -- | token parser for whitespace whiteSpace :: CharParser st () whiteSpace = PT.whiteSpace lexer ++-- | token parser for colon+colon:: CharParser st String+colon = PT.colon lexer+++ -- | token parser for semicolon semi :: CharParser st String semi = PT.semi lexer + -- | token parser for comma comma :: CharParser st String comma = PT.comma lexer +-- | token parser for comma separated list of items+commaSep :: CharParser st a -> CharParser st [a]+commaSep = PT.commaSep lexer++ -- | token parser for symbol symbol :: String -> CharParser st String symbol = PT.symbol lexer+++-- | token parser for symbol+identifier :: CharParser st String+identifier = PT.identifier lexer+++-- | token parser for symbol+operator:: CharParser st String+operator = PT.operator lexer+
src/simgi.hs view
@@ -24,7 +24,7 @@ -- imports import Prelude import System.IO-import System.Random+--import System.Random import System.Environment -- local imports@@ -43,9 +43,10 @@ main = -- process command line arguments- getArgs >>= process_commandline - >>= \((_,_),files) -> + getArgs >>= process_commandline initialModelState+ >>= \(state, files) -> + -- reject anything but a single input file if length files /= 1 then usage@@ -56,28 +57,29 @@ >>= \content -> -- parse input file- case runParser input_parser initialModelState "" content of+ case runParser input_parser state "" content of Left er -> putStrLn (show er) Right parsedState -> case check_input parsedState of Left err -> putStrLn err- Right _ -> newStdGen- >>= \gen -> + Right _ -> -- set up simuation let - rands = randomRs (0,1) gen :: [Double] initialOutput = create_initial_output parsedState initialState = create_initial_state parsedState - rands initialOutput+ initialOutput totalTime = maxTime parsedState - dataDumpIter = maxIter parsedState+ dataDumpIter = outputBufferSize parsedState outFile = outfileName parsedState in -- open output file openFile outFile WriteMode - >>= \handle -> startup_message+ >>= \handle -> ++ -- print initial startup info+ startup_message initialState -- ready to run the simulation >> gillespie_driver handle totalTime dataDumpIter
test/Makefile view
@@ -3,13 +3,37 @@ default: check -.PHONY: check clean-check:- make -C rpnstack_test- make -C reversible_test+.PHONY: check check_event_parser check_rpnstack check_reversible \+ check_reaction_parser check_irreversible clean -clean:- make -C rpnstack_test clean- make -C reversible_test clean+check_rpnstack:+ make -C $@+++check_reversible:+ make -C $@++check_irreversible:+ make -C $@ ++check_event_parser:+ make -C $@+++check_reaction_parser:+ make -C $@++check: check_rpnstack check_event_parser check_reaction_parser \+ check_reversible check_irreversible ++++clean: + make -C check_rpnstack clean+ make -C check_reversible clean+ make -C check_irreversible clean+ make -C check_event_parser clean+ make -C check_reaction_parser clean+
+ test/check_event_parser/EventParserTest.hs view
@@ -0,0 +1,248 @@+{-----------------------------------------------------------------+ + (c) 2009 Markus Dittrich + + This program is free software; you can redistribute it + and/or modify it under the terms of the GNU General Public + License Version 3 as published by the Free Software Foundation. + + This program is distributed in the hope that it will be useful,+ but WITHOUT ANY WARRANTY; without even the implied warranty of+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the+ GNU General Public License Version 3 for more details.+ + You should have received a copy of the GNU General Public + License along with this program; if not, write to the Free + Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA.++--------------------------------------------------------------------}++-- | this routine tests some aspects of the event block parsing+-- routines+module Main where+++-- imports +import Control.Monad.Writer+import qualified Data.Map as M+import Prelude+import System.Exit++-- local imports+import Engine+import GenericModel+import PrettyPrint+import InputParser+import TestHelpers+import TokenParser+++-----------------------------------------------------------------+-- +-- test data specifications+--+-----------------------------------------------------------------++-- | expected values for parser events. The first+-- entry in the tuple is of type Bool giving the expected+-- value of the evaluated trigger expression. The second+-- is a Map of expected molecule counts after the actions+-- have been evaluated+type ExpectedOutput = (Bool, M.Map String Int)++++-- | set up the test data+-- Format : (parse expression, expected result)+type TestCase = (String, ExpectedOutput)++++-- | simple tests without access to local variables+--- NOTE: for now we simply test if parsing succeeds+simpleEventParseTests :: [TestCase]+simpleEventParseTests = + [ ("{ z == 100 } => { x = x }", + (True, M.fromList [("x",1000),("y",2000),("z",100)]))++ , ("{ z== 100} => { x = x+y; y = y+1 }", + (True, M.fromList [("x",3000),("y",2001),("z",100)]))++ , ("{z==121 } => { y=10 }", + (False, M.fromList [("x",1000),("y",10),("z",100)]))++ , ("{z==101} => { x = x/2.0 }", + (False, M.fromList [("x",500),("y",2000),("z",100)]))++ , ("{ z == 100}=>{ z= 121; x = 500; y=2000;}", + (True, M.fromList [("x",500),("y",2000),("z",121)]))++ , ("{ y >= 121 }=>{z=1050}", + (True, M.fromList [("x",1000),("y",2000),("z",1050)]))++ , ("{ z <= 100} => { x= sqrt(4)*100; y = x^2/10}", + (True, M.fromList [("x",200),("y",4000),("z",100)]))++ , ("{ x >= 100} => { z= z; z = 2*z; z = 3*z;}", + (True, M.fromList [("x",1000),("y",2000),("z",600)]))++ , ("{ z < 100} => { x= 10; y = z}", + (False, M.fromList [("x",10),("y",100),("z",100)]))++ , ("{ x > 100} => { y= x+z}", + (True, M.fromList [("x",1000),("y",1100),("z",100)]))++ , ("{ y == 2000 } => { x=1001 }", + (True, M.fromList [("x",1001),("y",2000),("z",100)]))++ , ("{ z <= 302 } => { x = 10; y = 10; z= 200; z=10 }", + (True, M.fromList [("x",10),("y",10),("z",10)]))++ , ("{ y == 100 } => { x = 10; y = 20; }", + (False, M.fromList [("x",10),("y",20),("z",100)]))++ , ("{ z == 100 } => { x = 10; y = 20 }", + (True, M.fromList [("x",10),("y",20),("z",100)]))++ , ("{ z == 100 } => { x = 10; y = 20 + z }", + (True, M.fromList [("x",10),("y",120),("z",100)]))++ , ("{ x > 100 } => { x = 10*x+y; y = 20+z^2 }", + (True, M.fromList [("x",12000),("y",10020),("z",100)]))++ , ("{ z == 100 } => { x = 1e6*exp(-TIME); y = x; }",+ (True, M.fromList [("x",4),("y",4),("z",100)]))++ , ("{ y == 100 } => { x = 2; x = y^x; y = x*exp(-TIME) }", + (False, M.fromList [("x",4000000),("y",17),("z",100)]))++ , ("{ z == 100 } => { x = 10+z; y = sqrt(log((x^2))); }",+ (True, M.fromList [("x",110),("y",3),("z",100)])) + ]+++----------------------------------------------------------------+-- Molecule maps and definitions for specific simulation times+----------------------------------------------------------------++-- | testmap containing a set of molecules and their +-- concentrations+testMap_1 :: MoleculeMap+testMap_1 = M.fromList [("x",1000),("y",2000),("z",100)]++-- | a simulation time+time_1 :: Double+time_1 = 12.345+++++----------------------------------------------------------------+--+-- main driver routines+--+----------------------------------------------------------------+++-- | main test driver+main :: IO ()+main = putStrLn "\n\n\nTesting Input Parser"++ -- check event parser+ >> (putStr $ color_string Cyan "\nEvent parse tests:\n")+ >> let eventParseOut = execWriter $ event_parser_test_driver + testMap_1 time_1 simpleEventParseTests+ in+ examine_output eventParseOut >>= \eventParseStatus ->+++ -- check event triggers and actions+ (putStr $ color_string Cyan "\n\nEvent trigger/action tests:\n")+ >> let eventActionOut = execWriter $ event_action_test_driver + testMap_1 time_1 simpleEventParseTests+ in+ examine_output eventActionOut >>= \eventActionStatus ->+++ -- evaluate status and return+ let status = eventParseStatus && eventActionStatus + in+ if status == True then+ exitWith ExitSuccess+ else+ exitWith $ ExitFailure 1 ++++-- | this driver parses the event expression and checks that+-- there are no errors (there shouldn't be any)+event_parser_test_driver :: MoleculeMap -> Double -> [TestCase] + -> Writer [TestResult] ()+event_parser_test_driver _ _ [] = return ()+event_parser_test_driver mol t (x:xs) =++ let expr = fst x+ expected = snd x+ in++ -- parse expression+ case runParser parse_events testModelState "" expr of+ Left er -> tell [TestResult False expr (show expected) (show er)]+ Right _ -> + tell [TestResult True expr (show expected) ("good parse")]+ >> event_parser_test_driver mol t xs+ +++-- | this driver parses the event expression and checks if+-- the trigger and its actions action check out, i.e., we apply +-- them to our current state and see if we get the proper number +-- of molecules.+-- NOTE: We do not thread the state through all tests but always+-- start with our default, otherwise it becomes very tedious+-- to add/remove tests.+-- This parser should be run after the proper parsing has been+-- verified via event_parser_test_driver+event_action_test_driver :: MoleculeMap -> Double -> [TestCase] + -> Writer [TestResult] ()+event_action_test_driver _ _ [] = return ()+event_action_test_driver mol t (x:xs) =++ let expr = fst x+ expectedTrigger = fst . snd $ x+ expectedMols = snd . snd $ x + in++ -- parse expression+ case runParser parse_events testModelState "" expr of+ Left er -> tell [TestResult False expr "" (show er)]+ Right event -> + + -- make sure the trigger expression evaluated properly+ let+ actions = evtActions event+ outMols = execute_actions actions (SymbolTable mol M.empty) t + outTrigger = check_trigger mol t expectedTrigger event + in+ case outTrigger && ((molSymbols outMols) == expectedMols) of++ False -> tell [TestResult False expr + (show expectedTrigger ++ " => " ++ show expectedMols)+ (show outTrigger ++ " => " ++ show (molSymbols outMols))]++ True -> tell [TestResult True expr "" ("good parse")]+ >> event_action_test_driver mol t xs++++-- | given an Event and a MoleculeMap check that the trigger+-- evaluates to the expected value+check_trigger :: MoleculeMap -> Double -> Bool -> Event+ -> Bool+check_trigger molMap t expected (Event { evtTrigger = trigger }) =+ computed == expected+ + where+ computed = compute_trigger (SymbolTable molMap M.empty) t trigger+
+ test/check_event_parser/Makefile view
@@ -0,0 +1,15 @@+include ../test_helpers/make_common++default: check+++.PHONY: check clean++check:+ ghc $(GHC_FLAGS_DEVEL) -i$(SRC_LOCATION) -i$(HELPER_LOCATION) \+ --make EventParserTest.hs+ @./EventParserTest+++clean:+ rm -f *.hi *.o EventParserTest
+ test/check_irreversible/Makefile view
@@ -0,0 +1,19 @@+include ../test_helpers/make_common+++default: check++.PHONY: check clean+check:+ @echo+ @echo "Running irreversible reaction test 1"+ @echo "This may take a while ...."+ @./run_test.1.sh $(SIMGI_PATH)+ @echo "Running irreversible reaction test 2"+ @echo "This may take a while ...."+ @./run_test.2.sh $(SIMGI_PATH)+ @echo+++clean:+ rm -f *.dat
+ test/check_irreversible/irreversible.1.sgl view
@@ -0,0 +1,22 @@++def parameters+ time = 1.0+ systemVol = 1.25e-19 + outputBuffer = 50000 + outputFreq = 100+ outputFile = "irreversible.1.dat"+end++def molecules+ a = 1000+ b = 1000+ c = 0+end++def reactions+ a + b -> c | 1e6 |+end ++def output+ [a, b, c]+end
+ test/check_irreversible/irreversible.2.sgl view
@@ -0,0 +1,22 @@++def parameters+ time = 1.0+ systemVol = 1.25e-19 + outputBuffer = 50000 + outputFreq = 100+ outputFile = "irreversible.2.dat"+end++def molecules+ a = 1000+ b = 1000+ c = 0+end++def reactions+ b -> c | 1e6 |+end ++def output+ [a, b, c]+end
+ test/check_irreversible/run_test.1.sh view
@@ -0,0 +1,60 @@+#!/bin/bash+#+# short wrapper for runing a number of reversible+# simulations and then checking if we can reproduce+# the correct average product concentration+#++# expected values of molecule counts+a_expect=0.000000000000000+b_expect=0.000000000000000+c_expect=1000.000000000000000++# +simgi_exe="${1}"++# our global status+status=0++for ((counter=0; counter <= 100; counter++)); do++ # provide a little "progressbar"+ printf "."++ # run and process+ ${simgi_exe} irreversible.1.sgl >& /dev/null ++ a=$(tail -n 1 irreversible.1.dat | gawk ' { print $3 }')+ b=$(tail -n 1 irreversible.1.dat | gawk ' { print $4 }')+ c=$(tail -n 1 irreversible.1.dat | gawk ' { print $5 }')++ if [[ ${a_expect} != ${a} || ${b_expect} != ${b} \+ || ${c_expect} != ${c} ]];+ then+ status=1+ fi++ # unlink+ rm -f irreversible.1.dat || return 1+ +done++# check if all the tests passed++# brief output+if [[ ${status} == 0 ]]; then+ echo+ echo+ echo "Congratulations - the irreversible reaction test 1 passed!"+ echo+ echo+else+ echo+ echo+ echo "Error - the irreversible reaction test failed. Please check!"+ echo+ echo+fi++# return status of deviation check+exit ${status}
+ test/check_irreversible/run_test.2.sh view
@@ -0,0 +1,60 @@+#!/bin/bash+#+# short wrapper for runing a number of reversible+# simulations and then checking if we can reproduce+# the correct average product concentration+#++# expected values of molecule counts+a_expect=1000.000000000000000+b_expect=0.000000000000000+c_expect=1000.000000000000000++# +simgi_exe="${1}"++# our global status+status=0++for ((counter=0; counter <= 100; counter++)); do++ # provide a little "progressbar"+ printf "."++ # run and process+ ${simgi_exe} irreversible.2.sgl >& /dev/null ++ a=$(tail -n 1 irreversible.2.dat | gawk ' { print $3 }')+ b=$(tail -n 1 irreversible.2.dat | gawk ' { print $4 }')+ c=$(tail -n 1 irreversible.2.dat | gawk ' { print $5 }')++ if [[ ${a_expect} != ${a} || ${b_expect} != ${b} \+ || ${c_expect} != ${c} ]];+ then+ status=1+ fi++ # unlink+ rm -f irreversible.2.dat || return 1+ +done++# check if all the tests passed++# brief output+if [[ ${status} == 0 ]]; then+ echo+ echo+ echo "Congratulations - the irreversible reaction test 2 passed!"+ echo+ echo+else+ echo+ echo+ echo "Error - the irreversible reaction test failed. Please check!"+ echo+ echo+fi++# return status of deviation check+exit ${status}
+ test/check_reaction_parser/Makefile view
@@ -0,0 +1,15 @@+include ../test_helpers/make_common++default: check+++.PHONY: check clean++check:+ ghc $(GHC_FLAGS_DEVEL) -i$(SRC_LOCATION) -i$(HELPER_LOCATION) \+ --make ReactionParserTest.hs+ @./ReactionParserTest+++clean:+ rm -f *.hi *.o ReactionParserTest
+ test/check_reaction_parser/ReactionParserTest.hs view
@@ -0,0 +1,200 @@+{-----------------------------------------------------------------+ + (c) 2009 Markus Dittrich + + This program is free software; you can redistribute it + and/or modify it under the terms of the GNU General Public + License Version 3 as published by the Free Software Foundation. + + This program is distributed in the hope that it will be useful,+ but WITHOUT ANY WARRANTY; without even the implied warranty of+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the+ GNU General Public License Version 3 for more details.+ + You should have received a copy of the GNU General Public + License along with this program; if not, write to the Free + Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA.++--------------------------------------------------------------------}++-- | this routine tests some aspects of the reaction block+-- parsing routines+module Main where+++-- imports +import Control.Monad.Writer+import qualified Data.Map as M+import Prelude+import System.Exit++-- local imports+import GenericModel+import PrettyPrint+import InputParser+import RpnData+import TestHelpers+import TokenParser+++-----------------------------------------------------------------+-- +-- test data specifications+--+-----------------------------------------------------------------++-- | expected values for parser events. The first+-- entry in the tuple is of type Bool giving the expected+-- value of the evaluated trigger expression. The second+-- is a Map of expected molecule counts after the actions+-- have been evaluated+type ExpectedOutput = Reaction++++-- | set up the test data+-- Format : (parse expression, expected result)+type TestCase = (String, ExpectedOutput)++++-- | simple tests without access to local variables+--- NOTE: for now we simply test if parsing succeeds+simpleReactParseTests :: [TestCase]+simpleReactParseTests = + [ ("x + y -> z | 1e6 |",+ Reaction (Constant 1e6) + [("y",id),("x",id)]+ [("x",-1),("y",-1),("z",1)])++ , ("x + y + z -> w | 10.3 |",+ Reaction (Constant 10.3)+ [("z",id),("y",id),("x",id)]+ [("x",-1),("y",-1),("z",-1),("w",1)])++ , ("x + y + z -> w | {x*y} |",+ Reaction (Function . RpnStack $ [Variable "x", Variable "y", + BinFunc (*)])+ [("z",id),("y",id),("x",id)]+ [("x",-1),("y",-1),("z",-1),("w",1)])++ , ("2x + 2y + z -> nil | { exp(-TIME) } |",+ Reaction (Function . RpnStack $ [Variable "TIME", Number (-1.0)+ , BinFunc (*), UnaFunc exp])+ [("z",id),("y",\a -> 0.5 * a * (a-1))+ ,("x",\a -> 0.5 * a * (a-1))]+ [("x",-2),("y",-2),("z",-1)])++ , ("2x + y + z -> 2x + 20y + 5z | {x+y+z } |",+ Reaction (Function . RpnStack $ [Variable "x", Variable "y", + BinFunc (+), Variable "z", BinFunc (+)])+ [("z",id),("y",id),("x",\a -> 0.5 * a * (a-1))]+ [("y",19),("z",4)])++ , ("2x + y + z -> 2x + 20y + 5z | { x* y* z } |",+ Reaction (Function . RpnStack $ [Variable "x", Variable "y", + BinFunc (*), Variable "z", BinFunc (*)])+ [("z",id),("y",id),("x",\a -> 0.5 * a * (a-1))]+ [("y",19),("z",4)])++ , ("2 x+ y+ z-> 2 x + 20 y + 5z|{x*y * z } |",+ Reaction (Function . RpnStack $ [Variable "x", Variable "y", + BinFunc (*), Variable "z", BinFunc (*)])+ [("z",id),("y",id),("x",\a -> 0.5 * a * (a-1))]+ [("y",19),("z",4)])++ , ("x + y + z -> x + y + z |{ x-x+y-y+z-z }|",+ Reaction (Function . RpnStack $ [Number 0.0])+ [("z",id),("y",id),("x",id)]+ [])++ , ("nil -> x + y + z |{ x-x+y-y+z-z }|",+ Reaction (Function . RpnStack $ [Number 0.0])+ []+ [("x",1),("y",1),("z",1)])++ , ("nil -> nil |{ sqrt(2.0) } |",+ Reaction (Function . RpnStack $ [Number 2.0, UnaFunc sqrt])+ []+ [])+ ]++++----------------------------------------------------------------+-- Molecule maps and definitions for specific simulation times+----------------------------------------------------------------++-- | testmap containing a set of molecules and their +-- concentrations+testMap_1 :: MoleculeMap+testMap_1 = M.fromList [("x",1000),("y",2000),("z",100)]++-- | a simulation time+time_1 :: Double+time_1 = 12.345+++++----------------------------------------------------------------+--+-- main driver routines+--+----------------------------------------------------------------+++-- | main test driver+main :: IO ()+main = putStrLn "\n\n\nTesting Reaction Parser"++ -- check event parser+ >> (putStr $ color_string Cyan "\nReaction parse tests:\n")+ >> let reactParseOut = execWriter $ react_parser_test_driver + testMap_1 time_1 simpleReactParseTests+ in+ examine_output reactParseOut >>= \reactParseStatus ->+++ -- evaluate status and return+ let status = reactParseStatus + in+ if status == True then+ exitWith ExitSuccess+ else+ exitWith $ ExitFailure 1 ++++-- the trigger and its actions action check out, i.e., we apply +-- them to our current state and see if we get the proper number +-- of molecules.+-- NOTE: We do not thread the state through all tests but always+-- start with our default, otherwise it becomes very tedious+-- to add/remove tests.+-- This parser should be run after the proper parsing has been+-- verified via event_parser_test_driver+react_parser_test_driver :: MoleculeMap -> Double -> [TestCase] + -> Writer [TestResult] ()+react_parser_test_driver _ _ [] = return ()+react_parser_test_driver mol t (x:xs) =++ let+ parseString = fst x+ expected = snd x+ in++ -- parse expression+ case runParser parse_reaction testModelState "" parseString of+ Left er -> tell [TestResult False parseString "" (show er)]+ Right react -> + + case (react == expected) of++ False -> tell [TestResult False parseString "N/A" "N/A"]++ True -> tell [TestResult True parseString "" "good parse"]+ >> react_parser_test_driver mol t xs++
+ test/check_reversible/Makefile view
@@ -0,0 +1,19 @@+include ../test_helpers/make_common+++default: check++.PHONY: check clean+check:+ ghc $(GHC_FLAGS_DEVEL) -i$(SRC_LOCATION) --make check_deviation.hs+ ghc $(GHC_FLAGS_DEVEL) -i$(SRC_LOCATION) --make average.hs+ @echo+ @echo "Running reversible reaction test."+ @echo "This may take a while ...."+ @./run_test.sh $(SIMGI_PATH)+ @echo+ ++clean:+ rm -f *.hi *.o check_deviation average reversible.dat \+ simgi.*
+ test/check_reversible/average.hs view
@@ -0,0 +1,38 @@+{-----------------------------------------------------------------+ + (c) 2009 Markus Dittrich + + This program is free software; you can redistribute it + and/or modify it under the terms of the GNU General Public + License Version 3 as published by the Free Software Foundation. + + This program is distributed in the hope that it will be useful,+ but WITHOUT ANY WARRANTY; without even the implied warranty of+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the+ GNU General Public License Version 3 for more details.+ + You should have received a copy of the GNU General Public + License along with this program; if not, write to the Free + Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA.++--------------------------------------------------------------------}++-- | short helper script for averaging a file containing a single+-- column of doubles+module Main where+++-- imports+import Prelude+++main :: IO ()+main = getContents+ >>= \content -> + let+ items = map (read) . lines $ content :: [Double]+ theSum = foldr (+) 0 items+ average = theSum/(fromIntegral . length $ items)+ in+ putStrLn (show average)
+ test/check_reversible/check_deviation.hs view
@@ -0,0 +1,46 @@+{-----------------------------------------------------------------+ + (c) 2009 Markus Dittrich + + This program is free software; you can redistribute it + and/or modify it under the terms of the GNU General Public + License Version 3 as published by the Free Software Foundation. + + This program is distributed in the hope that it will be useful,+ but WITHOUT ANY WARRANTY; without even the implied warranty of+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the+ GNU General Public License Version 3 for more details.+ + You should have received a copy of the GNU General Public + License along with this program; if not, write to the Free + Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA.++--------------------------------------------------------------------}++-- | short helper script for averaging a file containing a single+-- column of doubles+module Main where++-- imports+import Prelude+import System+import System.Environment+++main :: IO Int+main = getArgs + >>= \(name:value:tol:_) -> readFile name+ >>= \content -> + let+ items = map (read) . lines $ content :: [Double]+ theSum = foldr (+) 0 items+ average = theSum/(fromIntegral . length $ items)+ expected = read value :: Double+ tolerance = read tol :: Double+ deviation = abs((average - expected)/expected)+ in+ if ( deviation < tolerance )+ then exitWith ExitSuccess+ else exitWith $ ExitFailure 1+
+ test/check_reversible/reversible.sgl view
@@ -0,0 +1,29 @@+{-----------------------------------------------------++ this is the input deck for the oregonator model + (C) 2009 Markus Dittrich++------------------------------------------------------}++def parameters+ time = 1.0e-4+ systemVol = 1.25e-19+ outputBuffer = 50000+ outputFreq = 10+ outputFile = "reversible.dat"+end++def molecules+ a = 1000+ b = 1000+ c = 0+end++def reactions+ a + b -> c | {1e6*10} |+ c -> a + b | 1e5 |+end++def output+ [a, b, c]+end
+ test/check_reversible/run_test.sh view
@@ -0,0 +1,53 @@+#!/bin/bash+#+# short wrapper for runing a number of reversible+# simulations and then checking if we can reproduce+# the correct average product concentration+#++# +simgi_exe="${1}"++# create a global tempfile to collect the data+globalFile=$(mktemp simgi.XXXXXXXXXXXXX)++for ((counter=0; counter <= 100; counter++)); do++ # provide a little "progressbar"+ printf "."++ # run and process+ ${simgi_exe} -s ${counter} reversible.sgl >& /dev/null + tail -n 400 reversible.dat | gawk ' { print $5}' | ./average >> ${globalFile} || return 1++ # unlink+ rm -f reversible.dat || return 1+ +done++# check if we're within the specs+# the expected number of products is 430.643462709951+# and we allow 0.5% tolerance+./check_deviation ${globalFile} 430.643462709951 5e-3+status=$?++# brief output+if [[ ${status} == 0 ]]; then+ echo+ echo+ echo "Congratulations - the reversible reaction test passed!"+ echo+ echo+else+ echo+ echo+ echo "Error - the reversible reaction test failed. Please check!"+ echo+ echo+fi++# remove files+rm -f ${globalFile}++# return status of deviation check+exit ${status}
+ test/check_rpnstack/Makefile view
@@ -0,0 +1,15 @@+include ../test_helpers/make_common+++default: check+++.PHONY: check clean+check:+ ghc $(GHC_FLAGS_DEVEL) -i$(SRC_LOCATION) -i$(HELPER_LOCATION) \+ --make RpnStackTest.hs+ @./RpnStackTest+++clean:+ rm -f *.hi *.o RpnStackTest
+ test/check_rpnstack/RpnStackTest.hs view
@@ -0,0 +1,182 @@+{-----------------------------------------------------------------+ + (c) 2009 Markus Dittrich + + This program is free software; you can redistribute it + and/or modify it under the terms of the GNU General Public + License Version 3 as published by the Free Software Foundation. + + This program is distributed in the hope that it will be useful,+ but WITHOUT ANY WARRANTY; without even the implied warranty of+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the+ GNU General Public License Version 3 for more details.+ + You should have received a copy of the GNU General Public + License along with this program; if not, write to the Free + Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA.++--------------------------------------------------------------------}++-- | this routine tests the functionality in our RPNStack, i.e.,+-- the infix to RPN parser as well as the compute engine+module Main where+++-- imports +import Control.Monad.Writer+import qualified Data.Map as M+import Prelude+import System.Exit++-- local imports+import ExtraFunctions+import GenericModel+import PrettyPrint+import RpnParser+import RpnCalc+import TestHelpers+import TokenParser+++-----------------------------------------------------------------+-- +-- test data specifications+--+-----------------------------------------------------------------++-- | set up the test data+-- Format : (parse expression, expected result)+type TestCase = (String, Double)++++-- | simple tests without access to local variables+simpleTests :: [TestCase]+simpleTests = + [ ("3+4", 7.0) + , ("3 +4", 7.0)+ , ("3+ 4", 7.0)+ , ("3 + 4 ", 7.0)+ , ("-3 + 3 ", 0.0)+ , ("3*4+4", 16.0)+ , ("3 * 4+4 ", 16.0)+ , ("3+4*5", 23.0)+ , ("(3+4)*5", 35.0)+ , ("( 3 + 4 )* 5", 35.0)+ , ("sqrt(2)^2", 2.0)+ , ("(-2)^2", 4.0)+ , ("exp(-1)", 0.36787944117144233)+ , ("log(exp(3)) * 3^2", 27.0)+ , ("log ( exp ( 3) ) *3 ^2", 27.0)+ , ("((((((((3*(3+(3*3) + 4) + 2) +3)-34) -4)+1))))", 16.0)+ , ("(1*(1 * (1 *(((((3*(3+(3*3) + 4) + 2) +3)-34) -4)+1))))", 16.0)+ , ("log(exp(2) - sqrt(2))", 1.7875577437560926)+ , ("2 * 3.14 * sqrt(2)", 8.88126117170303786)+ ]++++----------------------------------------------------------------+-- tests with access to local variables and time+----------------------------------------------------------------++-- | variable tests+variableTests :: [TestCase]+variableTests =+ [ ("3*x", 3000)+ , ("x", 1000)+ , ("sqrt(x)^2", 1000)+ , ("x+y*z", 1000)+ , ("x + y * z", 1000)+ , ("(x+y)*z", 0)+ , ("(x + y)* z", 0)+ , ("exp(z)*TIME", 12.345)+ , ("exp(z ) * TIME ", 12.345)+ , ("-x * -y", 2.0e6)+ , ("x*exp(-TIME)", 4.351456244655325e-3)+ , ("x-x +x -x -y + y", 0.0)+ , ("-TIME/TIME + TIME - TIME + 1.0", 0.0)+ ]++++----------------------------------------------------------------+-- Molecule maps and definitions for specific simulation times+----------------------------------------------------------------++-- | testmap containing a set of molecules and their +-- concentrations+testMap_1 :: MoleculeMap+testMap_1 = M.fromList [("x",1000),("y",2000),("z",0)]++-- | a simulation time+time_1 :: Double+time_1 = 12.345+++++----------------------------------------------------------------+--+-- main driver routines+--+----------------------------------------------------------------+++-- | main test driver+main :: IO ()+main = putStrLn "\n\n\nTesting RPN stack (parser/compute engine)"++ -- run simple tests+ >> (putStr $ color_string Cyan "\nSimple tests:\n")+ >> let simpleOut = execWriter $ test_driver + testMap_1 time_1 simpleTests+ in+ examine_output simpleOut >>= \simpleStatus ->+++ -- run variable tests+ (putStr $ color_string Cyan "\n\nVariable tests:\n")+ >> let varOut = execWriter $ test_driver testMap_1 time_1+ variableTests+ in+ examine_output varOut >>= \varStatus ->+++ -- evaluate status and return+ let status = simpleStatus && varStatus in+ if status == True then+ exitWith ExitSuccess+ else+ exitWith $ ExitFailure 1 +++-- | driver for running a test routine that results in a+-- successful evaluation of a test expression+test_driver :: MoleculeMap -> Double -> [TestCase] + -> Writer [TestResult] ()+test_driver _ _ [] = return ()+test_driver molMap t (x:xs) =++ let expr = fst x+ expected = snd x+ in++ -- parse expression+ case runParser parse_infix_to_rpn testModelState "" expr of+ Left er -> tell [TestResult False expr (show expected) (show er)]+ Right stack ->++ -- evalute RPN stack+ let result = rpn_compute (SymbolTable molMap M.empty) t stack in+ examine_result expected result expr + >> test_driver molMap t xs++ where+ examine_result target out anExpr = + if is_equal target out+ then + tell [TestResult True anExpr (show target) (show out)]+ else + tell [TestResult False anExpr (show target) (show out)]
− test/reversible_test/Makefile
@@ -1,19 +0,0 @@--SRC_LOCATION=../../src-SIMGI_PATH=${SRC_LOCATION}/simgi--default: check--.PHONY: check clean-check:- ghc -i$(SRC_LOCATION) --make check_deviation.hs- ghc -i$(SRC_LOCATION) --make average.hs- @echo- @echo "Running reversible reaction test."- @echo "This may take a while ...."- @./run_test.sh ${SIMGI_PATH}- @echo- --clean:- rm -f *.hi *.o check_deviation average
− test/reversible_test/average.hs
@@ -1,34 +0,0 @@-{------------------------------------------------------------------ - (c) 2009 Markus Dittrich - - This program is free software; you can redistribute it - and/or modify it under the terms of the GNU General Public - License Version 3 as published by the Free Software Foundation. - - This program is distributed in the hope that it will be useful,- but WITHOUT ANY WARRANTY; without even the implied warranty of- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the- GNU General Public License Version 3 for more details.- - You should have received a copy of the GNU General Public - License along with this program; if not, write to the Free - Software Foundation, Inc., 59 Temple Place - Suite 330, - Boston, MA 02111-1307, USA.----------------------------------------------------------------------}---- | short helper script for averaging a file containing a single--- column of doubles-module Main where---main :: IO ()-main = getContents- >>= \content -> - let- items = map (read) . lines $ content :: [Double]- theSum = foldr (+) 0 items- average = theSum/(fromIntegral . length $ items)- in- putStrLn (show average)
− test/reversible_test/check_deviation.hs
@@ -1,45 +0,0 @@-{------------------------------------------------------------------ - (c) 2009 Markus Dittrich - - This program is free software; you can redistribute it - and/or modify it under the terms of the GNU General Public - License Version 3 as published by the Free Software Foundation. - - This program is distributed in the hope that it will be useful,- but WITHOUT ANY WARRANTY; without even the implied warranty of- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the- GNU General Public License Version 3 for more details.- - You should have received a copy of the GNU General Public - License along with this program; if not, write to the Free - Software Foundation, Inc., 59 Temple Place - Suite 330, - Boston, MA 02111-1307, USA.----------------------------------------------------------------------}---- | short helper script for averaging a file containing a single--- column of doubles-module Main where---- imports-import System-import System.Environment---main :: IO Int-main = getArgs - >>= \(name:value:tol:_) -> readFile name- >>= \content -> - let- items = map (read) . lines $ content :: [Double]- theSum = foldr (+) 0 items- average = theSum/(fromIntegral . length $ items)- expected = read value :: Double- tolerance = read tol :: Double- deviation = abs((average - expected)/expected)- in- if ( deviation < tolerance )- then exitWith ExitSuccess- else exitWith $ ExitFailure 1-
− test/reversible_test/reversible.sgl
@@ -1,25 +0,0 @@-{------------------------------------------------------- this is the input deck for the oregonator model - (C) 2009 Markus Dittrich--------------------------------------------------------}--def parameters- time = 1.0e-4- systemVol = 1.25e-19- outputIter = 50000- outputFreq = 10- outputFile = "reversible.dat"-end--def molecules- a 1000- b 1000- c 0-end--def reactions- a + b -> c { 1e6*10 }- c -> a + b { 1e5 }-end
− test/reversible_test/run_test.sh
@@ -1,53 +0,0 @@-#!/bin/bash-#-# short wrapper for runing a number of reversible-# simulations and then checking if we can reproduce-# the correct average product concentration-#--# -simgi_exe="${1}"--# create a global tempfile to collect the data-globalFile=$(mktemp simgi.XXXXXXXXXXXXX)--for ((counter=0; counter <= 100; counter++)); do-- # provide a little "progressbar"- printf "."-- # run and process- ${simgi_exe} reversible.sgl >& /dev/null - tail -n 400 reversible.dat | gawk ' { print $5}' | ./average >> ${globalFile} || return 1-- # unlink- rm -f reversible.dat || return 1- -done--# check if we're within the specs-# the expected number of products is 430.643462709951-# and we allow 0.5% tolerance-./check_deviation ${globalFile} 430.643462709951 5e-3-status=$?--# brief output-if [[ ${status} == 0 ]]; then- echo- echo- echo "Congratulations - the reversible reaction test passed!"- echo- echo-else- echo- echo- echo "Error - the reversible reaction test failed. Please check!"- echo- echo-fi--# remove files-rm -f ${globalFile}--# return status of deviation check-exit ${status}
− test/rpnstack_test/Makefile
@@ -1,15 +0,0 @@--SRC_LOCATION=../../src---default: check---.PHONY: check clean-check:- ghc -i$(SRC_LOCATION) --make RpnStackTest.hs- @./RpnStackTest---clean:- rm -f *.hi *.o RpnStackTest
− test/rpnstack_test/RpnStackTest.hs
@@ -1,227 +0,0 @@-{------------------------------------------------------------------ - (c) 2009 Markus Dittrich - - This program is free software; you can redistribute it - and/or modify it under the terms of the GNU General Public - License Version 3 as published by the Free Software Foundation. - - This program is distributed in the hope that it will be useful,- but WITHOUT ANY WARRANTY; without even the implied warranty of- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the- GNU General Public License Version 3 for more details.- - You should have received a copy of the GNU General Public - License along with this program; if not, write to the Free - Software Foundation, Inc., 59 Temple Place - Suite 330, - Boston, MA 02111-1307, USA.----------------------------------------------------------------------}---- | this routine tests the functionality in our RPNStack, i.e.,--- the infix to RPN parser as well as the compute engine-module Main where----- imports -import Control.Monad.Writer-import qualified Data.Map as M-import Prelude-import System.Exit---- local imports-import ExtraFunctions-import GenericModel-import PrettyPrint-import RpnParser-import RpnCalc-import TokenParser----------------------------------------------------------------------- --- test data specifications------------------------------------------------------------------------- | set up the test data--- Format : (parse expression, expected result)-type TestCase = (String, Double)------ | simple tests without access to local variables-simpleTests :: [TestCase]-simpleTests = - [ ("3+4", 7.0) - , ("3 +4", 7.0)- , (" 3+ 4", 7.0)- , (" 3 + 4 ", 7.0)- , ("-3 + 3 ", 0.0)- , ("3*4+4", 16.0)- , ("3 * 4+4 ", 16.0)- , ("3+4*5", 23.0)- , ("(3+4)*5", 35.0)- , (" ( 3 + 4 )* 5", 35.0)- , ("sqrt(2)^2", 2.0)- , ("(-2)^2", 4.0)- , ("exp(-1)", 0.36787944117144233)- , ("log(exp(3)) * 3^2", 27.0)- , ("log ( exp ( 3) ) *3 ^2", 27.0)- , ("((((((((3*(3+(3*3) + 4) + 2) +3)-34) -4)+1))))", 16.0)- , ("(1*(1 * (1 *(((((3*(3+(3*3) + 4) + 2) +3)-34) -4)+1))))", 16.0)- , ("log(exp(2) - sqrt(2))", 1.7875577437560926)- , ("2 * 3.14 * sqrt(2)", 8.88126117170303786)- ]----- | variable tests-variableTests :: [TestCase]-variableTests =- [ ("3*x", 3000)- , ("sqrt(x)^2", 1000)- , ("x+y*z", 1000)- , ("x + y * z", 1000)- , ("(x+y)*z", 0)- , ("(x + y)* z", 0)- , ("exp(z)*TIME", 12.345)- , ("exp(z ) * TIME ", 12.345)- , ("-x * -y", 2.0e6)- , ("x*exp(-TIME)", 4.351456244655325e-3)- , ("x-x +x -x -y + y", 0.0)- , ("-TIME/TIME + TIME - TIME + 1.0", 0.0)- ]---------------------------------------------------------------------- tests with access to local variables and time---------------------------------------------------------------------------------------------------------------------------------------- Molecule maps and definitions for specific simulation times--------------------------------------------------------------------- | testmap containing a set of molecules and their --- concentrations-testMap_1 :: MoleculeMap-testMap_1 = M.fromList [("x",1000),("y",2000),("z",0)]---- | a simulation time-time_1 :: Double-time_1 = 12.345--------------------------------------------------------------------------- main driver routines------------------------------------------------------------------------- | main test driver-main :: IO ()-main = putStrLn "\n\n\nTesting RPN stack (parser/compute engine)"-- -- run simple tests- >> (putStr $ color_string Cyan "\nSimple tests:\n")- >> let simpleOut = execWriter $ test_driver - testMap_1 time_1 simpleTests- in- examine_output simpleOut >>= \simpleStatus ->--- -- run variable tests- (putStr $ color_string Cyan "\n\nVariable tests:\n")- >> let varOut = execWriter $ test_driver testMap_1 time_1- variableTests- in- examine_output varOut >>= \varStatus ->--- -- evaluate status and return- let status = simpleStatus && varStatus in- if status == True then- exitWith ExitSuccess- else- exitWith $ ExitFailure 1 ----- | examine the output of a test routine--- | helper function for examining the output of a good test run--- (i.e. one that should succeed), prints out the result for each --- test, collects the number of successes/failures and returns --- True in case all tests succeeded and False otherwise-examine_output :: [TestResult] -> IO Bool-examine_output = foldM examine_output_h True- - where- examine_output_h :: Bool -> TestResult -> IO Bool- examine_output_h acc (TestResult status token target actual) = do- if status == True then do- putStr $ color_string Blue "["- putStr $ color_string White "OK"- putStr $ color_string Blue "] "- putStr $ color_string Green " Successfully evaluated "- putStrLn $ color_string Yellow token- return $ acc && True- else do- putStr $ color_string Blue "["- putStr $ color_string Red "TROUBLE"- putStr $ color_string Blue "] "- putStr $ color_string Green " Failed to evaluate "- putStrLn $ color_string Yellow token- putStrLn $ color_string Green "\t\texpected : " - ++ (show target)- putStrLn $ color_string Green "\t\tgot : " - ++ (show actual)- return False----- | driver for running a test routine that results in a--- successful evaluation of a test expression-test_driver :: MoleculeMap -> Double -> [TestCase] - -> Writer [TestResult] ()-test_driver _ _ [] = return ()-test_driver mol time (x:xs) =-- let expr = fst x- expected = snd x- in-- -- parse expression- case runParser parse_infix_to_rpn initialModelState "" expr of- Left er -> tell [TestResult False expr (show expected) (show er)]- Right stack ->-- -- evalute RPN stack- let result = rpn_compute mol time stack in- examine_result expected result expr - >> test_driver mol time xs-- where- examine_result target out expr = - if is_equal target out- then - tell [TestResult True expr (show target) (show out)]- else - tell [TestResult False expr (show target) (show out)]- ---- | data structure for keeping track of our test results--- which consist of a bool indicating success--- or failure, the test token as well as the expected and--- received result-data TestResult = TestResult { status :: Bool- , token :: String- , target :: String- , actual :: String- }---defaultResult :: TestResult-defaultResult = TestResult False "" "" ""-
+ test/test_helpers/TestHelpers.hs view
@@ -0,0 +1,110 @@+{-----------------------------------------------------------------+ + (c) 2009 Markus Dittrich + + This program is free software; you can redistribute it + and/or modify it under the terms of the GNU General Public + License Version 3 as published by the Free Software Foundation. + + This program is distributed in the hope that it will be useful,+ but WITHOUT ANY WARRANTY; without even the implied warranty of+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the+ GNU General Public License Version 3 for more details.+ + You should have received a copy of the GNU General Public + License along with this program; if not, write to the Free + Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA.++--------------------------------------------------------------------}++-- | this module provides a few common routines used in our+-- unit tests+module TestHelpers ( examine_output+ , testModelState+ , TestResult(TestResult)+ , defaultResult+ ) where+++-- imports +import Control.Monad+import qualified Data.Map as M+import Prelude++-- local imports+import GenericModel+import PrettyPrint++++-- | examine the output of a test routine+-- | helper function for examining the output of a good test run+-- (i.e. one that should succeed), prints out the result for each +-- test, collects the number of successes/failures and returns +-- True in case all tests succeeded and False otherwise+examine_output :: [TestResult] -> IO Bool+examine_output = foldM examine_output_h True+ + where+ examine_output_h :: Bool -> TestResult -> IO Bool+ examine_output_h acc (TestResult state tok targ act) = do+ if state == True then do+ putStr $ color_string Blue "["+ putStr $ color_string White "OK"+ putStr $ color_string Blue "] "+ putStr $ color_string Green " Successfully evaluated "+ putStrLn $ color_string Yellow tok+ return $ acc && True+ else do+ putStr $ color_string Blue "["+ putStr $ color_string Red "TROUBLE"+ putStr $ color_string Blue "] "+ putStr $ color_string Green " Failed to evaluate "+ putStrLn $ color_string Yellow tok+ putStrLn $ color_string Green "\t\texpected : " + ++ (show targ)+ putStrLn $ color_string Green "\t\tgot : " + ++ (show act)+ return False++++-- | data structure for keeping track of our test results+-- which consist of a bool indicating success+-- or failure, the test token as well as the expected and+-- received result+data TestResult = TestResult { status :: Bool+ , token :: String+ , target :: String+ , actual :: String+ }+++defaultResult :: TestResult+defaultResult = TestResult { status = False + , token = ""+ , target = ""+ , actual = ""+ }+++-- | initial model state we use for our unit tests+-- we use a negative system volume so we don't have+-- to deal with unit conversion when testing the+-- reaction parser+testModelState :: ModelState+testModelState = ModelState { molCount = M.empty+ , rates = []+ , reactions = []+ , randNums = []+ , events = []+ , systemVol = -1.0 -- negative!+ , currentTime = 0.0+ , currentIter = 0+ , maxTime = 0.0+ , outputBufferSize = 10000+ , outputFreq = 1000+ , outputCache = []+ , outfileName = ""+ }
+ test/test_helpers/make_common view
@@ -0,0 +1,7 @@+# some common variables all Makefiles in the test directory use++GHC_FLAGS_DEVEL = -O -Wall -fwarn-simple-patterns -fwarn-tabs -fwarn-incomplete-record-updates -fwarn-monomorphism-restriction -fwarn-implicit-prelude -fno-warn-orphans++SRC_LOCATION=../../src+HELPER_LOCATION=../test_helpers+SIMGI_PATH=$(SRC_LOCATION)/simgi