diff --git a/.gitignore b/.gitignore new file mode 100644 --- /dev/null +++ b/.gitignore @@ -0,0 +1,11 @@ +syntax: glob +src/*.o +src/*.hi +src/simgi +src/*.swp +test/*.o +test/*.hi +test/*.swp +average +RpnStackTest +check_deviation diff --git a/AUTHORS b/AUTHORS new file mode 100644 --- /dev/null +++ b/AUTHORS @@ -0,0 +1,1 @@ +Markus Dittrich diff --git a/COPYING b/COPYING new file mode 100644 --- /dev/null +++ b/COPYING @@ -0,0 +1,674 @@ + GNU GENERAL PUBLIC LICENSE + Version 3, 29 June 2007 + + Copyright (C) 2007 Free Software Foundation, Inc. + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The GNU General Public License is a free, copyleft license for +software and other kinds of works. + + The licenses for most software and other practical works are designed +to take away your freedom to share and change the works. By contrast, +the GNU General Public License is intended to guarantee your freedom to +share and change all versions of a program--to make sure it remains free +software for all its users. We, the Free Software Foundation, use the +GNU General Public License for most of our software; it applies also to +any other work released this way by its authors. You can apply it to +your programs, too. + + When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +them if you wish), that you receive source code or can get it if you +want it, that you can change the software or use pieces of it in new +free programs, and that you know you can do these things. + + To protect your rights, we need to prevent others from denying you +these rights or asking you to surrender the rights. Therefore, you have +certain responsibilities if you distribute copies of the software, or if +you modify it: responsibilities to respect the freedom of others. + + For example, if you distribute copies of such a program, whether +gratis or for a fee, you must pass on to the recipients the same +freedoms that you received. You must make sure that they, too, receive +or can get the source code. And you must show them these terms so they +know their rights. + + Developers that use the GNU GPL protect your rights with two steps: +(1) assert copyright on the software, and (2) offer you this License +giving you legal permission to copy, distribute and/or modify it. + + For the developers' and authors' protection, the GPL clearly explains +that there is no warranty for this free software. For both users' and +authors' sake, the GPL requires that modified versions be marked as +changed, so that their problems will not be attributed erroneously to +authors of previous versions. + + Some devices are designed to deny users access to install or run +modified versions of the software inside them, although the manufacturer +can do so. This is fundamentally incompatible with the aim of +protecting users' freedom to change the software. The systematic +pattern of such abuse occurs in the area of products for individuals to +use, which is precisely where it is most unacceptable. Therefore, we +have designed this version of the GPL to prohibit the practice for those +products. If such problems arise substantially in other domains, we +stand ready to extend this provision to those domains in future versions +of the GPL, as needed to protect the freedom of users. + + Finally, every program is threatened constantly by software patents. +States should not allow patents to restrict development and use of +software on general-purpose computers, but in those that do, we wish to +avoid the special danger that patents applied to a free program could +make it effectively proprietary. To prevent this, the GPL assures that +patents cannot be used to render the program non-free. + + The precise terms and conditions for copying, distribution and +modification follow. + + TERMS AND CONDITIONS + + 0. Definitions. + + "This License" refers to version 3 of the GNU General Public License. + + "Copyright" also means copyright-like laws that apply to other kinds of +works, such as semiconductor masks. + + "The Program" refers to any copyrightable work licensed under this +License. Each licensee is addressed as "you". "Licensees" and +"recipients" may be individuals or organizations. + + To "modify" a work means to copy from or adapt all or part of the work +in a fashion requiring copyright permission, other than the making of an +exact copy. The resulting work is called a "modified version" of the +earlier work or a work "based on" the earlier work. + + A "covered work" means either the unmodified Program or a work based +on the Program. + + To "propagate" a work means to do anything with it that, without +permission, would make you directly or secondarily liable for +infringement under applicable copyright law, except executing it on a +computer or modifying a private copy. Propagation includes copying, +distribution (with or without modification), making available to the +public, and in some countries other activities as well. + + To "convey" a work means any kind of propagation that enables other +parties to make or receive copies. Mere interaction with a user through +a computer network, with no transfer of a copy, is not conveying. + + An interactive user interface displays "Appropriate Legal Notices" +to the extent that it includes a convenient and prominently visible +feature that (1) displays an appropriate copyright notice, and (2) +tells the user that there is no warranty for the work (except to the +extent that warranties are provided), that licensees may convey the +work under this License, and how to view a copy of this License. If +the interface presents a list of user commands or options, such as a +menu, a prominent item in the list meets this criterion. + + 1. Source Code. + + The "source code" for a work means the preferred form of the work +for making modifications to it. "Object code" means any non-source +form of a work. + + A "Standard Interface" means an interface that either is an official +standard defined by a recognized standards body, or, in the case of +interfaces specified for a particular programming language, one that +is widely used among developers working in that language. + + The "System Libraries" of an executable work include anything, other +than the work as a whole, that (a) is included in the normal form of +packaging a Major Component, but which is not part of that Major +Component, and (b) serves only to enable use of the work with that +Major Component, or to implement a Standard Interface for which an +implementation is available to the public in source code form. A +"Major Component", in this context, means a major essential component +(kernel, window system, and so on) of the specific operating system +(if any) on which the executable work runs, or a compiler used to +produce the work, or an object code interpreter used to run it. + + The "Corresponding Source" for a work in object code form means all +the source code needed to generate, install, and (for an executable +work) run the object code and to modify the work, including scripts to +control those activities. However, it does not include the work's +System Libraries, or general-purpose tools or generally available free +programs which are used unmodified in performing those activities but +which are not part of the work. For example, Corresponding Source +includes interface definition files associated with source files for +the work, and the source code for shared libraries and dynamically +linked subprograms that the work is specifically designed to require, +such as by intimate data communication or control flow between those +subprograms and other parts of the work. + + The Corresponding Source need not include anything that users +can regenerate automatically from other parts of the Corresponding +Source. + + The Corresponding Source for a work in source code form is that +same work. + + 2. Basic Permissions. + + All rights granted under this License are granted for the term of +copyright on the Program, and are irrevocable provided the stated +conditions are met. This License explicitly affirms your unlimited +permission to run the unmodified Program. The output from running a +covered work is covered by this License only if the output, given its +content, constitutes a covered work. This License acknowledges your +rights of fair use or other equivalent, as provided by copyright law. + + You may make, run and propagate covered works that you do not +convey, without conditions so long as your license otherwise remains +in force. You may convey covered works to others for the sole purpose +of having them make modifications exclusively for you, or provide you +with facilities for running those works, provided that you comply with +the terms of this License in conveying all material for which you do +not control copyright. Those thus making or running the covered works +for you must do so exclusively on your behalf, under your direction +and control, on terms that prohibit them from making any copies of +your copyrighted material outside their relationship with you. + + Conveying under any other circumstances is permitted solely under +the conditions stated below. Sublicensing is not allowed; section 10 +makes it unnecessary. + + 3. Protecting Users' Legal Rights From Anti-Circumvention Law. + + No covered work shall be deemed part of an effective technological +measure under any applicable law fulfilling obligations under article +11 of the WIPO copyright treaty adopted on 20 December 1996, or +similar laws prohibiting or restricting circumvention of such +measures. + + When you convey a covered work, you waive any legal power to forbid +circumvention of technological measures to the extent such circumvention +is effected by exercising rights under this License with respect to +the covered work, and you disclaim any intention to limit operation or +modification of the work as a means of enforcing, against the work's +users, your or third parties' legal rights to forbid circumvention of +technological measures. + + 4. Conveying Verbatim Copies. + + You may convey verbatim copies of the Program's source code as you +receive it, in any medium, provided that you conspicuously and +appropriately publish on each copy an appropriate copyright notice; +keep intact all notices stating that this License and any +non-permissive terms added in accord with section 7 apply to the code; +keep intact all notices of the absence of any warranty; and give all +recipients a copy of this License along with the Program. + + You may charge any price or no price for each copy that you convey, +and you may offer support or warranty protection for a fee. + + 5. Conveying Modified Source Versions. + + You may convey a work based on the Program, or the modifications to +produce it from the Program, in the form of source code under the +terms of section 4, provided that you also meet all of these conditions: + + a) The work must carry prominent notices stating that you modified + it, and giving a relevant date. + + b) The work must carry prominent notices stating that it is + released under this License and any conditions added under section + 7. This requirement modifies the requirement in section 4 to + "keep intact all notices". + + c) You must license the entire work, as a whole, under this + License to anyone who comes into possession of a copy. This + License will therefore apply, along with any applicable section 7 + additional terms, to the whole of the work, and all its parts, + regardless of how they are packaged. This License gives no + permission to license the work in any other way, but it does not + invalidate such permission if you have separately received it. + + d) If the work has interactive user interfaces, each must display + Appropriate Legal Notices; however, if the Program has interactive + interfaces that do not display Appropriate Legal Notices, your + work need not make them do so. + + A compilation of a covered work with other separate and independent +works, which are not by their nature extensions of the covered work, +and which are not combined with it such as to form a larger program, +in or on a volume of a storage or distribution medium, is called an +"aggregate" if the compilation and its resulting copyright are not +used to limit the access or legal rights of the compilation's users +beyond what the individual works permit. Inclusion of a covered work +in an aggregate does not cause this License to apply to the other +parts of the aggregate. + + 6. Conveying Non-Source Forms. + + You may convey a covered work in object code form under the terms +of sections 4 and 5, provided that you also convey the +machine-readable Corresponding Source under the terms of this License, +in one of these ways: + + a) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by the + Corresponding Source fixed on a durable physical medium + customarily used for software interchange. + + b) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by a + written offer, valid for at least three years and valid for as + long as you offer spare parts or customer support for that product + model, to give anyone who possesses the object code either (1) a + copy of the Corresponding Source for all the software in the + product that is covered by this License, on a durable physical + medium customarily used for software interchange, for a price no + more than your reasonable cost of physically performing this + conveying of source, or (2) access to copy the + Corresponding Source from a network server at no charge. + + c) Convey individual copies of the object code with a copy of the + written offer to provide the Corresponding Source. This + alternative is allowed only occasionally and noncommercially, and + only if you received the object code with such an offer, in accord + with subsection 6b. + + d) Convey the object code by offering access from a designated + place (gratis or for a charge), and offer equivalent access to the + Corresponding Source in the same way through the same place at no + further charge. You need not require recipients to copy the + Corresponding Source along with the object code. If the place to + copy the object code is a network server, the Corresponding Source + may be on a different server (operated by you or a third party) + that supports equivalent copying facilities, provided you maintain + clear directions next to the object code saying where to find the + Corresponding Source. Regardless of what server hosts the + Corresponding Source, you remain obligated to ensure that it is + available for as long as needed to satisfy these requirements. + + e) Convey the object code using peer-to-peer transmission, provided + you inform other peers where the object code and Corresponding + Source of the work are being offered to the general public at no + charge under subsection 6d. + + A separable portion of the object code, whose source code is excluded +from the Corresponding Source as a System Library, need not be +included in conveying the object code work. + + A "User Product" is either (1) a "consumer product", which means any +tangible personal property which is normally used for personal, family, +or household purposes, or (2) anything designed or sold for incorporation +into a dwelling. In determining whether a product is a consumer product, +doubtful cases shall be resolved in favor of coverage. For a particular +product received by a particular user, "normally used" refers to a +typical or common use of that class of product, regardless of the status +of the particular user or of the way in which the particular user +actually uses, or expects or is expected to use, the product. A product +is a consumer product regardless of whether the product has substantial +commercial, industrial or non-consumer uses, unless such uses represent +the only significant mode of use of the product. + + "Installation Information" for a User Product means any methods, +procedures, authorization keys, or other information required to install +and execute modified versions of a covered work in that User Product from +a modified version of its Corresponding Source. The information must +suffice to ensure that the continued functioning of the modified object +code is in no case prevented or interfered with solely because +modification has been made. + + If you convey an object code work under this section in, or with, or +specifically for use in, a User Product, and the conveying occurs as +part of a transaction in which the right of possession and use of the +User Product is transferred to the recipient in perpetuity or for a +fixed term (regardless of how the transaction is characterized), the +Corresponding Source conveyed under this section must be accompanied +by the Installation Information. But this requirement does not apply +if neither you nor any third party retains the ability to install +modified object code on the User Product (for example, the work has +been installed in ROM). + + The requirement to provide Installation Information does not include a +requirement to continue to provide support service, warranty, or updates +for a work that has been modified or installed by the recipient, or for +the User Product in which it has been modified or installed. Access to a +network may be denied when the modification itself materially and +adversely affects the operation of the network or violates the rules and +protocols for communication across the network. + + Corresponding Source conveyed, and Installation Information provided, +in accord with this section must be in a format that is publicly +documented (and with an implementation available to the public in +source code form), and must require no special password or key for +unpacking, reading or copying. + + 7. Additional Terms. + + "Additional permissions" are terms that supplement the terms of this +License by making exceptions from one or more of its conditions. +Additional permissions that are applicable to the entire Program shall +be treated as though they were included in this License, to the extent +that they are valid under applicable law. If additional permissions +apply only to part of the Program, that part may be used separately +under those permissions, but the entire Program remains governed by +this License without regard to the additional permissions. + + When you convey a copy of a covered work, you may at your option +remove any additional permissions from that copy, or from any part of +it. (Additional permissions may be written to require their own +removal in certain cases when you modify the work.) You may place +additional permissions on material, added by you to a covered work, +for which you have or can give appropriate copyright permission. + + Notwithstanding any other provision of this License, for material you +add to a covered work, you may (if authorized by the copyright holders of +that material) supplement the terms of this License with terms: + + a) Disclaiming warranty or limiting liability differently from the + terms of sections 15 and 16 of this License; or + + b) Requiring preservation of specified reasonable legal notices or + author attributions in that material or in the Appropriate Legal + Notices displayed by works containing it; or + + c) Prohibiting misrepresentation of the origin of that material, or + requiring that modified versions of such material be marked in + reasonable ways as different from the original version; or + + d) Limiting the use for publicity purposes of names of licensors or + authors of the material; or + + e) Declining to grant rights under trademark law for use of some + trade names, trademarks, or service marks; or + + f) Requiring indemnification of licensors and authors of that + material by anyone who conveys the material (or modified versions of + it) with contractual assumptions of liability to the recipient, for + any liability that these contractual assumptions directly impose on + those licensors and authors. + + All other non-permissive additional terms are considered "further +restrictions" within the meaning of section 10. If the Program as you +received it, or any part of it, contains a notice stating that it is +governed by this License along with a term that is a further +restriction, you may remove that term. If a license document contains +a further restriction but permits relicensing or conveying under this +License, you may add to a covered work material governed by the terms +of that license document, provided that the further restriction does +not survive such relicensing or conveying. + + If you add terms to a covered work in accord with this section, you +must place, in the relevant source files, a statement of the +additional terms that apply to those files, or a notice indicating +where to find the applicable terms. + + Additional terms, permissive or non-permissive, may be stated in the +form of a separately written license, or stated as exceptions; +the above requirements apply either way. + + 8. Termination. + + You may not propagate or modify a covered work except as expressly +provided under this License. Any attempt otherwise to propagate or +modify it is void, and will automatically terminate your rights under +this License (including any patent licenses granted under the third +paragraph of section 11). + + However, if you cease all violation of this License, then your +license from a particular copyright holder is reinstated (a) +provisionally, unless and until the copyright holder explicitly and +finally terminates your license, and (b) permanently, if the copyright +holder fails to notify you of the violation by some reasonable means +prior to 60 days after the cessation. + + Moreover, your license from a particular copyright holder is +reinstated permanently if the copyright holder notifies you of the +violation by some reasonable means, this is the first time you have +received notice of violation of this License (for any work) from that +copyright holder, and you cure the violation prior to 30 days after +your receipt of the notice. + + Termination of your rights under this section does not terminate the +licenses of parties who have received copies or rights from you under +this License. If your rights have been terminated and not permanently +reinstated, you do not qualify to receive new licenses for the same +material under section 10. + + 9. Acceptance Not Required for Having Copies. + + You are not required to accept this License in order to receive or +run a copy of the Program. Ancillary propagation of a covered work +occurring solely as a consequence of using peer-to-peer transmission +to receive a copy likewise does not require acceptance. However, +nothing other than this License grants you permission to propagate or +modify any covered work. These actions infringe copyright if you do +not accept this License. Therefore, by modifying or propagating a +covered work, you indicate your acceptance of this License to do so. + + 10. Automatic Licensing of Downstream Recipients. + + Each time you convey a covered work, the recipient automatically +receives a license from the original licensors, to run, modify and +propagate that work, subject to this License. You are not responsible +for enforcing compliance by third parties with this License. + + An "entity transaction" is a transaction transferring control of an +organization, or substantially all assets of one, or subdividing an +organization, or merging organizations. If propagation of a covered +work results from an entity transaction, each party to that +transaction who receives a copy of the work also receives whatever +licenses to the work the party's predecessor in interest had or could +give under the previous paragraph, plus a right to possession of the +Corresponding Source of the work from the predecessor in interest, if +the predecessor has it or can get it with reasonable efforts. + + You may not impose any further restrictions on the exercise of the +rights granted or affirmed under this License. For example, you may +not impose a license fee, royalty, or other charge for exercise of +rights granted under this License, and you may not initiate litigation +(including a cross-claim or counterclaim in a lawsuit) alleging that +any patent claim is infringed by making, using, selling, offering for +sale, or importing the Program or any portion of it. + + 11. Patents. + + A "contributor" is a copyright holder who authorizes use under this +License of the Program or a work on which the Program is based. The +work thus licensed is called the contributor's "contributor version". + + A contributor's "essential patent claims" are all patent claims +owned or controlled by the contributor, whether already acquired or +hereafter acquired, that would be infringed by some manner, permitted +by this License, of making, using, or selling its contributor version, +but do not include claims that would be infringed only as a +consequence of further modification of the contributor version. For +purposes of this definition, "control" includes the right to grant +patent sublicenses in a manner consistent with the requirements of +this License. + + Each contributor grants you a non-exclusive, worldwide, royalty-free +patent license under the contributor's essential patent claims, to +make, use, sell, offer for sale, import and otherwise run, modify and +propagate the contents of its contributor version. + + In the following three paragraphs, a "patent license" is any express +agreement or commitment, however denominated, not to enforce a patent +(such as an express permission to practice a patent or covenant not to +sue for patent infringement). To "grant" such a patent license to a +party means to make such an agreement or commitment not to enforce a +patent against the party. + + If you convey a covered work, knowingly relying on a patent license, +and the Corresponding Source of the work is not available for anyone +to copy, free of charge and under the terms of this License, through a +publicly available network server or other readily accessible means, +then you must either (1) cause the Corresponding Source to be so +available, or (2) arrange to deprive yourself of the benefit of the +patent license for this particular work, or (3) arrange, in a manner +consistent with the requirements of this License, to extend the patent +license to downstream recipients. "Knowingly relying" means you have +actual knowledge that, but for the patent license, your conveying the +covered work in a country, or your recipient's use of the covered work +in a country, would infringe one or more identifiable patents in that +country that you have reason to believe are valid. + + If, pursuant to or in connection with a single transaction or +arrangement, you convey, or propagate by procuring conveyance of, a +covered work, and grant a patent license to some of the parties +receiving the covered work authorizing them to use, propagate, modify +or convey a specific copy of the covered work, then the patent license +you grant is automatically extended to all recipients of the covered +work and works based on it. + + A patent license is "discriminatory" if it does not include within +the scope of its coverage, prohibits the exercise of, or is +conditioned on the non-exercise of one or more of the rights that are +specifically granted under this License. You may not convey a covered +work if you are a party to an arrangement with a third party that is +in the business of distributing software, under which you make payment +to the third party based on the extent of your activity of conveying +the work, and under which the third party grants, to any of the +parties who would receive the covered work from you, a discriminatory +patent license (a) in connection with copies of the covered work +conveyed by you (or copies made from those copies), or (b) primarily +for and in connection with specific products or compilations that +contain the covered work, unless you entered into that arrangement, +or that patent license was granted, prior to 28 March 2007. + + Nothing in this License shall be construed as excluding or limiting +any implied license or other defenses to infringement that may +otherwise be available to you under applicable patent law. + + 12. No Surrender of Others' Freedom. + + If conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot convey a +covered work so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you may +not convey it at all. For example, if you agree to terms that obligate you +to collect a royalty for further conveying from those to whom you convey +the Program, the only way you could satisfy both those terms and this +License would be to refrain entirely from conveying the Program. + + 13. Use with the GNU Affero General Public License. + + Notwithstanding any other provision of this License, you have +permission to link or combine any covered work with a work licensed +under version 3 of the GNU Affero General Public License into a single +combined work, and to convey the resulting work. The terms of this +License will continue to apply to the part which is the covered work, +but the special requirements of the GNU Affero General Public License, +section 13, concerning interaction through a network will apply to the +combination as such. + + 14. Revised Versions of this License. + + The Free Software Foundation may publish revised and/or new versions of +the GNU General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + + Each version is given a distinguishing version number. If the +Program specifies that a certain numbered version of the GNU General +Public License "or any later version" applies to it, you have the +option of following the terms and conditions either of that numbered +version or of any later version published by the Free Software +Foundation. If the Program does not specify a version number of the +GNU General Public License, you may choose any version ever published +by the Free Software Foundation. + + If the Program specifies that a proxy can decide which future +versions of the GNU General Public License can be used, that proxy's +public statement of acceptance of a version permanently authorizes you +to choose that version for the Program. + + Later license versions may give you additional or different +permissions. However, no additional obligations are imposed on any +author or copyright holder as a result of your choosing to follow a +later version. + + 15. Disclaimer of Warranty. + + THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY +APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT +HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY +OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, +THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM +IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF +ALL NECESSARY SERVICING, REPAIR OR CORRECTION. + + 16. Limitation of Liability. + + IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS +THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY +GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE +USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF +DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD +PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), +EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF +SUCH DAMAGES. + + 17. Interpretation of Sections 15 and 16. + + If the disclaimer of warranty and limitation of liability provided +above cannot be given local legal effect according to their terms, +reviewing courts shall apply local law that most closely approximates +an absolute waiver of all civil liability in connection with the +Program, unless a warranty or assumption of liability accompanies a +copy of the Program in return for a fee. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +state the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + + Copyright (C) + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + +Also add information on how to contact you by electronic and paper mail. + + If the program does terminal interaction, make it output a short +notice like this when it starts in an interactive mode: + + Copyright (C) + This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type `show c' for details. + +The hypothetical commands `show w' and `show c' should show the appropriate +parts of the General Public License. Of course, your program's commands +might be different; for a GUI interface, you would use an "about box". + + You should also get your employer (if you work as a programmer) or school, +if any, to sign a "copyright disclaimer" for the program, if necessary. +For more information on this, and how to apply and follow the GNU GPL, see +. + + The GNU General Public License does not permit incorporating your program +into proprietary programs. If your program is a subroutine library, you +may consider it more useful to permit linking proprietary applications with +the library. If this is what you want to do, use the GNU Lesser General +Public License instead of this License. But first, please read +. diff --git a/ChangeLog b/ChangeLog new file mode 100644 --- /dev/null +++ b/ChangeLog @@ -0,0 +1,8 @@ +2009-06-02 Markus Dittrich + + * 0.1.1 release (identical to 0.1 with only a few minor + changes needed for uploading to Hackage). + +2009-05-30 Markus Dittrich + + * initial commit diff --git a/INSTALL b/INSTALL new file mode 100644 --- /dev/null +++ b/INSTALL @@ -0,0 +1,16 @@ +COMPILATION +----------- +Compilaton of simgi requires + +- >=ghc-6.10 +- >=gmp-4.3 + +To compile the documentation (not required), you will also need + +- >=docutils-0.5 +- latex, e.g., tetex or texlive + +Building of simgi can be done either via + +- the standard ``make, make check, make install`` +- or via cabal diff --git a/Makefile b/Makefile new file mode 100644 --- /dev/null +++ b/Makefile @@ -0,0 +1,55 @@ +# Copyright 2009 Markus Dittrich +# Distributed under the terms of the GNU General Public License v3 + +VERSION=0.1.1 +DESTDIR= +prefix=/usr +mandir=$(DESTDIR)$(prefix)/share/man/man1 +docdir=$(DESTDIR)$(prefix)/share/doc/simgi-$(VERSION) +htmldir=$(docdir)/html +bindir=$(DESTDIR)$(prefix)/bin + +GHC_FLAGS_DEVEL = -O -Wall -fwarn-simple-patterns -fwarn-tabs -fwarn-incomplete-record-updates -fwarn-monomorphism-restriction -fwarn-implicit-prelude -fno-warn-orphans +GHC_FLAGS_RELEASE = -O2 + +OBJECTS = src/simgi.hs src/CommandLine.hs src/Engine.hs \ + src/GenericModel.hs src/InputCheck.hs src/InputParser.hs \ + src/PrettyPrint.hs src/RpnParser.hs src/RpnCalc.hs \ + src/RpnData.hs src/TokenParser.hs + +all: simgi + +simgi: $(OBJECTS) + ghc -i./Models -i./src $(GHC_FLAGS_RELEASE) --make src/simgi.hs + + + +debug: $(OBJECTS) + ghc -i./Models -i./src $(GHC_FLAGS_DEVEL) --make src/simgi.hs + + + +check: $(OBJECTS) all + 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/*.html $(htmldir)/ + + +.PHONY: clean doc + +doc: + make -C doc + +clean: + rm -f src/*.o src/*.hi src/simgi + make -C test clean + make -C doc clean diff --git a/Models/brusselator.sgl b/Models/brusselator.sgl new file mode 100644 --- /dev/null +++ b/Models/brusselator.sgl @@ -0,0 +1,26 @@ +{----------------------------------------------------- + + this is the input deck for the brusselator model + (C) 2009 Markus Dittrich + +------------------------------------------------------} + +def parameters + time = 50.0 + outputIter = 50000 + outputFreq = 200 + systemVol = nil -- interpret rates as propensities + outputFile = "brusselator_output.dat" +end + +def molecules + x 1000 + y 2000 +end + +def reactions + nil -> x { 5000 } + x -> y { 50.0 } + 2x + y -> 3x { 0.00005 } + x -> nil { 5.0 } +end diff --git a/Models/oregonator.sgl b/Models/oregonator.sgl new file mode 100644 --- /dev/null +++ b/Models/oregonator.sgl @@ -0,0 +1,28 @@ +{----------------------------------------------------- + + this is the input deck for the oregonator model + (C) 2009 Markus Dittrich + +------------------------------------------------------} + +def parameters + time = 50.0 + outputIter = 50000 + outputFreq = 500 + systemVol = nil -- interpret rates as propensities + outputFile = "oregonator_output.dat" +end + +def molecules + 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 diff --git a/Models/volterra.sgl b/Models/volterra.sgl new file mode 100644 --- /dev/null +++ b/Models/volterra.sgl @@ -0,0 +1,26 @@ +{----------------------------------------------------- + + this is the input deck for the lotka-volterra model + (C) 2009 Markus Dittrich + +------------------------------------------------------} + +def parameters + time = 50.0 + outputIter = 50000 + outputFreq = 200 + systemVol = nil -- interpret rates as propensities + outputFile = "volterra_output.dat" +end + +def molecules + x 1000 + y 2000 +end + +def reactions + x -> 2x { 10.0 } + x + y -> 2y { 0.01 } + y -> nil { 10.0 } +end + diff --git a/README b/README new file mode 100644 --- /dev/null +++ b/README @@ -0,0 +1,13 @@ +DESCRIPTION: +------------ + +simgi is a stochastic simulator using the Gillespie algorithm [1,2]. +simgi is released under the GPL v3. + +(C) 2009 Markus Dittrich + +See doc/simgi.html for more information. + + +[1] D.T. Gillespie, J. Comput. Phys., 22, 403 (1976) +[2] D.T. Gillespie, J. Phys. Chem., 81, 2340 (1977) diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,3 @@ +#!/usr/bin/env runhaskell +import Distribution.Simple +main = defaultMain diff --git a/doc/Makefile b/doc/Makefile new file mode 100644 --- /dev/null +++ b/doc/Makefile @@ -0,0 +1,31 @@ + +RST_FILES = simgi.rst + +default: doc + + +doc: clean_old build_html build_pdf + + +build_html: $(RST_FILES) + rst2html.py $< > $(<:.rst=.html) + +build_pdf: $(RST_FILES) + rst2latex.py $< > $(<:.rst=.tex) + + # fix indentation + sed -i "17 s:^\(.*\):\\\setlength{\\\parindent}{0in}\n\1:" \ + $(<:.rst=.tex) + pdflatex $(<:.rst=.tex) + rm -f *.aux *.log *.out *.tex + + +clean_old: + rm -f *.html *.pdf + + +.PHONY: clean + +clean: + + diff --git a/doc/simgi.html b/doc/simgi.html new file mode 100644 --- /dev/null +++ b/doc/simgi.html @@ -0,0 +1,535 @@ + + + + + + +simgi - A Stochastic Gillespie Simulator for Molecular Systems + + + + +
+

simgi - A Stochastic Gillespie Simulator for Molecular Systems

+ +++ + + + + + + + +
Author:Markus Dittrich
email:haskelladdict at users dot sourceforge dot net
Version:0.1.1 (06/02/2009)
+ +
+

Introduction

+

simgi is a fairly simple and straightforward stochastic simulator +based on Gillspie'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.

+
+
+

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.

+
+
+

Download

+

The current release of simgi can be downloaded here.

+
+
+

Compilation

+

Compilaton of simgi requires

+ +

To compile the documentation (not required), you will also need

+ +

Building of simgi can be done either via

+
    +
  • the standard make, make check, make install
    • +
  • or via cabal
    • +
+
+
+

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.

+

A SGL file consists of zero or more descriptor blocks of the form

+
+def <block name>
+
+  <block content>
+
+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. +However, it is strongly recommended to stick to a consistent and +"visually simple" layout to aid in "comprehending" the underlying +model.

+

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 +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 +respective block commands and block content:

+

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>
+
Maximum simulation time in seconds. Default is 0.0 s.
+
outputIter = <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. +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 +simulation (see outputFreq parameter).

+
+
outputFreq = <Integer>
+
Frequency with which output is generated. Default is 1000.
+
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>:
+
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

+
+

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
+end
+
+
+

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 }
+
+

Here, reactants and products are of the form

+
+<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 +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.

+

Here 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) }
+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 +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.

+
+
+
+

Example Input Files

+

Below are several example input files detailing the use of SGL:

+ +

These are also available in the Models/ sub-directory in the source tree.

+
+
+

Bugs

+

Please report all bugs and feature requests to +<haskelladdict at users dot sourceforge dot net>.

+
+
+

References

+ + + + + +
[1](1, 2) 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
+ + + + + +
[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.
+
+
+ + diff --git a/doc/simgi.pdf b/doc/simgi.pdf new file mode 100644 Binary files /dev/null and b/doc/simgi.pdf differ diff --git a/doc/simgi.rst b/doc/simgi.rst new file mode 100644 --- /dev/null +++ b/doc/simgi.rst @@ -0,0 +1,248 @@ +================================================================ +simgi - A Stochastic Gillespie Simulator for Molecular Systems +================================================================ + +:Author: Markus Dittrich + +:email: haskelladdict at users dot sourceforge dot net + +:Version: 0.1.1 (06/02/2009) + + +Contents +--------- + +1) Introduction_ +2) Status_ +3) Download_ +4) Compilation_ +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 +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 `_. + + +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. + + +Download +-------- + +The current release of simgi can be downloaded `here `_. + + +Compilation +----------- + +Compilaton of **simgi** requires + +- `>=ghc-6.10 `_ +- `>=gmp-4.3 `_ + +To compile the documentation (not required), you will also need + +- `>=docutils-0.5 `_ +- latex, e.g., tetex or texlive + + +Building of **simgi** can be done either via + +- the standard ``make, make check, make install`` +- or via cabal + + + +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. + +A SGL file consists of zero or more descriptor blocks of the form + +:: + + def + + + + 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. +However, it is strongly recommended to stick to a consistent and +"visually simple" layout to aid in "comprehending" the underlying +model. + +**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 + 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 +respective block commands and block content: + + +**parameter block:** * = parameters* + + The purpose of the parameter block is to describe the global + simulation parameters. The following parameters are currently + supported: + + *time = * + Maximum simulation time in seconds. Default is 0.0 s. + + *outputIter = * + 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. + 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 + simulation (see outputFreq parameter). + + *outputFreq = * + Frequency with which output is generated. Default is 1000. + + *systemVol = * + 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 = *: + Name of the output file. This is the only required parameter + in the parameter section. If not given, the simulation will + terminate. + + + +**molecule block:** * = molecules* + + This block consist of a list of pairs of the form :: + + + + 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 + end + + +**reaction block**: * = reactions* + + This block describes the reactions between molecules defined in + the molecule block. Reactions are specified via :: + + reactants -> product { rate expression } + + Here, ``reactants`` and ``products`` are of the form :: + + + + ..... + + In this expression, ```` is the reactant or product name + as defined in the molecule block and ```` an optional + integer specifying the stoichiometry. If ```` is not + explicitly given, it is assumed to be 1. + + The reaction rate can either be a fixed value of type ```` + or else an mathematical expression involving ````, + 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. + + Here 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) } + 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 + 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. + + +Example Input Files +------------------- + +Below are several example input files detailing the use of SGL: + +- `Lotka-Volterra Model `_ +- `Brusselator Model `_ +- `Oregonator Model `_ + +These are also available in the *Models/* sub-directory in the source tree. + + +Bugs +---- + +Please report all bugs and feature requests to +. + + +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 + +.. [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``. diff --git a/simgi.cabal b/simgi.cabal new file mode 100644 --- /dev/null +++ b/simgi.cabal @@ -0,0 +1,24 @@ +Name: simgi +Version: 0.1.1 +License: GPL +license-file: COPYING +copyright: (C) 2009 Markus Dittrich +category: Scientific Simulation +Synopsis: stochastic simulation engine +Description: simgi is a stochastic simulation engine to model + molecular systems using Gillespie's method. +Author: +Maintainer: +stability: alpha +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, + parsec == 2.1.*, mtl >= 1.1.0.2, haskell98, + random >= 1.0.0.1 + ghc-options: -O2 + Main-Is: simgi.hs + hs-source-dirs: src diff --git a/src/CommandLine.hs b/src/CommandLine.hs new file mode 100644 --- /dev/null +++ b/src/CommandLine.hs @@ -0,0 +1,85 @@ +{----------------------------------------------------------------- + + (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. + +--------------------------------------------------------------------} + +-- | main gsim driver +module CommandLine ( process_commandline + , Options(..) + ) where + +-- imports +import Prelude +import System +import System.Console.GetOpt + +-- local imports +import Messages + + +-- | main driver for command line processing +process_commandline :: [String] -> IO ((CmdlRequest,String), [String]) +process_commandline args = + + let + (actions, nonOpts, _) = getOpt RequireOrder options args + in + foldl (>>=) ( return defaultOptions ) actions >>= \opts -> + + let + Options { cmdlRequest = request + , cmdlString = pattern + } = opts + in + return ((request,pattern),nonOpts) + + + +-- | possible options for commandline +data CmdlRequest = None | Help + + +-- | data structure for keeping track of +-- selected command line options +data Options = Options { + cmdlRequest :: CmdlRequest, + cmdlString :: String +} + + +-- | default selections +defaultOptions :: Options +defaultOptions = Options { + cmdlRequest = None, + cmdlString = "" +} + + +-- | available command line flags +options :: [OptDescr (Options -> IO Options)] +options = [ + Option ['v'] ["version-info"] (NoArg version_info) + "show version information"] + + +-- | extractor function for version info +version_info :: Options -> IO Options +version_info _ = + do + show_version + exitWith ExitSuccess diff --git a/src/Engine.hs b/src/Engine.hs new file mode 100644 --- /dev/null +++ b/src/Engine.hs @@ -0,0 +1,210 @@ +{----------------------------------------------------------------- + + (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. + +--------------------------------------------------------------------} + +-- | the main compute Engine +module Engine ( create_initial_output + , create_initial_state + , gillespie_driver + , module GenericModel + ) where + +-- imports +import Control.Monad.State +import qualified Data.Map as M +import Prelude +import Text.Printf +import System.IO + +-- local imports +import ExtraFunctions +import GenericModel +import RpnCalc + +-- import Debug.Trace + +-- | main simulation driver +-- the simulator either stops when +-- 1) the number of iterations is exhausted +-- 2) the current time is > t_max, if t_max is set to +-- 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 + in + (write_output handle . reverse $ output) + >> if curTime >= simTime + then return () + else gillespie_driver handle simTime dmpIter newState + + + +-- | 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 = [] }) + + +-- | 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 + }) -> + + -- 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 + } + 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! + if (is_equal tau 0.0) + then let finalState = newState { currentTime = t_max } in + put finalState >> return output + else + if ( it_max == it || t >= t_max ) + then return output + else put newState >> run_gillespie + + +-- | 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 + | mod it afreq /= 0 = outlist + | otherwise = new_out:outlist + + where + new_out = Output { iteration = it + , time = t + , mols = amol + } + + + +-- | depending on which reaction happened adjust the number of +-- molecules in the system +adjust_molcount :: MoleculeMap -> [Reaction] -> Int -> MoleculeMap +adjust_molcount theMap rs mID = + + let (Reaction { react = 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 + m_new = M.insert k (a+val) m + in + adjustMap rands 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 ( MoleculeMap -> Double + -> RateList -> RateList +compute_rates [] _ _ rts = reverse rts +compute_rates ((Reaction {rate = c_in, aList = a_in }):rs) + theMap theTime rts = + + case c_in of + (Constant aRate) -> compute_rates rs theMap theTime + ((a_new aRate): rts) + (Function rateFunc) -> compute_rates rs theMap theTime + ((a_new . (rpn_compute theMap theTime) $ rateFunc):rts) + + where + mult = product $ map (\(a,f) -> f . fromIntegral $ + (M.!) theMap a) a_in + a_new = (*) mult + + +-- | initialize the output data structure +create_initial_output :: ModelState -> Output +create_initial_output (ModelState { molCount = initialMols }) = + + Output { iteration = 1 + , time = 0.0 + , mols = initialMols + } + + +-- | set up the initial state +create_initial_state:: ModelState -> [Double] -> Output -> ModelState +create_initial_state state rand output = + + state { rates = defaultRateList + , randNums = rand + , currentTime = 0.0 + , currentIter = 1 + , outputList = [output] + } + + +-- | 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 + in + -- write molecule data to output file + hPutStrLn handle (header ++ counts) + + -- write current iteration to stdout + >> (putStrLn $ "iteration | time ---> " ++ header) + >> write_output handle xs + + where + create_count_string :: MoleculeMap -> String + create_count_string = foldr (\x a -> (printf "%10d " x) ++ a) + "" . M.elems diff --git a/src/ExtraFunctions.hs b/src/ExtraFunctions.hs new file mode 100644 --- /dev/null +++ b/src/ExtraFunctions.hs @@ -0,0 +1,153 @@ +{-# LANGUAGE ForeignFunctionInterface #-} +{----------------------------------------------------------------- + + (c) 2008-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. + +--------------------------------------------------------------------} + +-- | definition of additional math and helper functions +module ExtraFunctions ( avogadroNum + , erf + , erfc + , fact + , is_equal + , is_equal_with + , real_exp + , to_int + , to_positive_int + ) where + + +-- imports +import Foreign() +import Foreign.C.Types +import Prelude + + +-- | a few constants +avogadroNum :: Double +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 +-- non-negative integers such as factorial. +-- To check if we are dealing with Double, we convert to an +-- 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 = + 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 = + 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., +-- pow a x = exp $ x * log a +foreign import ccall "math.h pow" + c_pow :: CDouble -> CDouble -> CDouble + +real_exp :: Double -> Double -> Double +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) +-- we use a recursive solution of the Taylor series expansion +erf :: Double -> Double +erf x + | x == 0.0 = 0.0 -- our recursive alg. loops forever + -- in this case + | abs(x) > 2.2 = 1.0 - erfc x -- use erfc for numerical accuracy + | otherwise = 2.0 / sqrt(pi) * (erf_h x x 1.0) + + where + erf_h old x_old n = let x_new = x_old * (x_next n) + tot = old + x_new + in + if abs(x_new/tot) < dbl_epsilon + then tot + else erf_h tot x_new (n+1.0) + + -- Note: We need (2::Int) here to silence ghc + x_next n = -(x^(2::Int)) * (2.0*n-1.0)/(n * (2.0*n+1.0)) + + + +-- | complementary error function erfc(x) = 1 - erf(x) +-- we use a recursive solution of the continued fraction +-- expression of erfc(x) for it superior convergence +-- property. Here, we calculate the ith and (i+1)th convergent, (see +-- http://archives.math.utk.edu/articles/atuyl/confrac/intro.html) +-- and terminate when the relative difference is smaller than a +-- certain threshold. +erfc :: Double -> Double +erfc x + | abs(x) < 2.2 = 1.0 - erf(x) -- use erf(x) in [-2.2,2.2] + | signum(x) < 0 = 2.0 - erfc(-x) -- continued fraction expansion + -- only valid for x > 0 + | otherwise = 1/sqrt(pi) * exp(-x^(2::Int)) + * (erfc_h nc1 nc2 dc1 dc2 1.0) + + where + nc1 = 1.0 :: Double -- numerator of 1st convergent + nc2 = x :: Double -- numerator of 2nd convergent + dc1 = x :: Double -- denominator of 1st convergent + dc2 = x^(2::Int)+0.5 :: Double -- denominator of 2nd convergent + + erfc_h n1 n2 d1 d2 i = + let num_new = n1*i + n2*x + denom_new = d1*i + d2*x + d_old = n2/d2 + d_new = num_new/denom_new + in + if abs((d_old - d_new)/d_new) < dbl_epsilon + then d_new + else erfc_h n2 num_new d2 denom_new (i+0.5) + diff --git a/src/GenericModel.hs b/src/GenericModel.hs new file mode 100644 --- /dev/null +++ b/src/GenericModel.hs @@ -0,0 +1,119 @@ +{----------------------------------------------------------------- + + (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. + +--------------------------------------------------------------------} + +-- | data structures needed for defining a stochastic model +module GenericModel ( defaultRateList + , GillespieState + , initialModelState + , ModelState(..) + , MoleculeMap + , Output(..) + , Rate(..) + , RateList + , Reaction(..) + ) where + + +-- imports +import Control.Monad.State +import qualified Data.Map as M +import Prelude + + +-- local imports +import RpnData + + +-- | 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 + +-- | List of reactions and corresponding rates +type RateList = [Double] + +defaultRateList :: RateList +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 +data Reaction = Reaction { rate :: Rate + , aList :: [(String,Double -> Double)] + , react :: [(String,Int)] + } + + +-- | 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 + } + +type GillespieState a = State ModelState a + + +-- | data structure for keeping track of our output +data Output = Output { iteration :: Integer + , time :: Double + , mols :: MoleculeMap + } + 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 = "" + } + + diff --git a/src/InputCheck.hs b/src/InputCheck.hs new file mode 100644 --- /dev/null +++ b/src/InputCheck.hs @@ -0,0 +1,128 @@ +{----------------------------------------------------------------- + + (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. + +--------------------------------------------------------------------} + +-- | module responsible for doing some basic input checking +module InputCheck ( check_input ) where + +-- imports +import Control.Monad.Error() +import qualified Data.List as L +import qualified Data.Map as M +import Prelude + +-- local imports +import GenericModel +import RpnData + + +-- | small function checking for common errors in the input deck +-- since checking is currently fairly simple we'll do it by hand +-- 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_molecules (defined_mols 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) + + where + -- | extract all reaction participants + react_mols = L.nub . L.concat . map (map (fst) . react) + + + -- | extract all definied molecules + defined_mols = M.keys + + + -- | 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 + + 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 + + +-- | make sure the user specified an output file name +check_filename :: String -> Either String Bool +check_filename name + | name == "" = Left "Error: Please specify an output file name!" + | otherwise = Right True + + +-- | make sure all molecules in reactions are defined +check_molecules :: [String] -> [String] -> Either String Bool +check_molecules defMols reactMols = + let + not_present = reactMols L.\\ defMols + in + case null not_present of + True -> Right True + False -> Left $ + "Error: The following molecules are not definied: " + ++ (L.concat $ L.intersperse "," not_present) + + +-- | make sure the user entered a positive integer for outputFreq +check_positive_outfreq :: Integer -> Either String Bool +check_positive_outfreq freq = + if freq < 0 + then Left "Error: outputFreq must be a positive integer!" + else Right True + + +-- | make sure the user entered a positive integer for outputFreq +check_positive_itercount :: Integer -> Either String Bool +check_positive_itercount iter = + if iter < 0 + then Left "Error: outputIter must be a positive integer!" + else Right True + + +-- | 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 = + let + no_mol = rateMols L.\\ defMols + in + case null no_mol of + True -> Right True + False -> Left $ + "Error: The following molecules defined in reaction " + ++ "rates do not exist: " + ++ (L.concat $ L.intersperse "," no_mol) diff --git a/src/InputParser.hs b/src/InputParser.hs new file mode 100644 --- /dev/null +++ b/src/InputParser.hs @@ -0,0 +1,317 @@ +{----------------------------------------------------------------- + + (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. + +--------------------------------------------------------------------} + +-- | input file parser +module InputParser ( input_parser ) where + +-- imports +import Control.Monad +import qualified Data.Map as M +import Prelude +import TokenParser + + +-- local imports +import ExtraFunctions +import GenericModel +import RpnData +import RpnParser + +-- import Debug.Trace + + +-- | main parser entry point +input_parser :: CharParser ModelState ModelState +input_parser = whiteSpace + *> many ( choice [ try parse_parameter_def + , try parse_molecule_def + , try parse_reaction_def + ]) + *> eof + >> getState + "main parser" + + +-- | parser for simulation parameters +parse_parameter_def :: CharParser ModelState () +parse_parameter_def = between (reserved "def" *> reserved "parameters") + (reserved "end") + (parse_parameters `sepBy` whiteSpace) + *> pure () + "parameter definitions" + + +-- | parse the individual parameters +parse_parameters :: CharParser ModelState () +parse_parameters = parse_time + <|> parse_outputFile + <|> parse_outputIter + <|> parse_outputFreq + <|> parse_systemVol + "time, outputIter, systemVol, outputFreq, outputFile" + + +-- | parse the simulation time specs +parse_time :: CharParser ModelState () +parse_time = join (updateState <$> insert_time + <$> (reserved "time" *> reservedOp "=" + *> parse_number)) + + where + 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 + <$> (reserved "systemVol" *> reservedOp "=" + *> (parse_positive_number + <|> parse_systemVol_nil ))) + "system volume" + + where + -- needed to avoid monomorphism warning + parse_systemVol_nil :: CharParser ModelState Double + parse_systemVol_nil = reserved "nil" *> pure (-1.0) + + 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 () +parse_outputFile = join (updateState <$> insert_filename + <$> (reserved "outputFile" *> reservedOp "=" + *> parse_filename )) + + where + insert_filename name state = state { outfileName = name } + + +-- | 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 "=" + *> integer )) + + where + insert_outputIter i state = state { maxIter = i } + + +-- | parse the output iteration specification if present +parse_outputFreq :: CharParser ModelState () +parse_outputFreq = join (updateState <$> insert_outputFreq + <$> (reserved "outputFreq" *> reservedOp "=" + *> integer )) + + where + 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) ) + "molecule definitions" + + where + insert_molecules :: [(String, Int)] -> ModelState -> ModelState + insert_molecules theMols state = + 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 + where + make_molecule mol aCount = (mol,fromInteger aCount) + + +-- | parser for a molecule name +-- 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 + + +-- | 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) ) + "reaction definitions" + + where + insert_reactions :: [Reaction] -> ModelState -> ModelState + 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 +-- by the user) we also need to extract the system +-- volume) +parse_reaction :: CharParser ModelState Reaction +parse_reaction = setup_reaction + <$> (parse_react_prod <* reservedOp "->") + <*> parse_react_prod + <*> parse_rate + <*> (getState + >>= \(ModelState {systemVol = vol}) -> pure vol) + where + -- | set up a Reaction data structure from the parsed reaction + setup_reaction r p cin vol = + let + action = create_react r p + hFactor = create_hFact r + theRate = if (vol < 0.0) -- no rate conversion for systemVol = nil + then cin + else convert_rate cin (M.size r) vol + in + Reaction { rate = theRate + , aList = hFactor + , react = action + } + + + -- | convert reaction propensities into rates if requested + -- by the user. For constants we simply multiply, for + -- rate functions we push the neccessary conversion onto + -- the stack + convert_rate theConst@(Constant c) order volume = + case order of + 1 -> theConst + _ -> Constant $ c/(avogadroNum * volume^(order-1)) + + convert_rate theFunc@(Function stack) order volume = + case order of + 1 -> theFunc + _ -> let mult = 1.0/(avogadroNum * volume^(order-1)) in + Function $ stack ++ [Number mult,BinFunc (*)] + + + + -- | create the list holding the molecule number changes for + -- this reaction + create_react r p = let + reacts = M.map (*(-1)) r + in + M.assocs $ 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 + -- 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 + + where + create_hFact_h acc [] = acc + create_hFact_h acc ((k,v):xs) = + let + v_int = fromIntegral v :: Double + in + create_hFact_h ((k,\x -> (1.0/v_int) + * generate_lambda v_int x):acc) xs + + where + generate_lambda :: Double -> Double -> Double + generate_lambda 1 x = x + 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)) + `sepBy` reservedOp "+") ) + "reactant or product list" + + where + 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 +parse_number = converter <$> naturalOrFloat + "signed integer or double" + where + converter val = case val of + Left i -> (fromInteger i) + 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 +parse_positive_number = naturalOrFloat + >>= \num -> case num of + Left ival -> if (ival > 0) + then return (fromInteger ival) + else pzero + Right dval -> if (dval > 0) + then return dval + else pzero + + "unsigned integer or double" diff --git a/src/Messages.hs b/src/Messages.hs new file mode 100644 --- /dev/null +++ b/src/Messages.hs @@ -0,0 +1,45 @@ +{----------------------------------------------------------------- + + (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. + +--------------------------------------------------------------------} + +-- | main gsim driver +module Messages ( show_version + , startup_message + , usage + ) where + +-- imports +import Prelude + + +-- | show version info +show_version :: IO () +show_version = putStrLn "This is simgi v0.1.1 (C) 2009 Markus Dittrich" + + +-- | show a brief startup message +startup_message :: IO () +startup_message = show_version + >> putStrLn "\nstarting simulation ..... here we go\n" + + + +-- | provide brief usage info +usage :: IO () +usage = putStrLn "Usage: simgi " diff --git a/src/PrettyPrint.hs b/src/PrettyPrint.hs new file mode 100644 --- /dev/null +++ b/src/PrettyPrint.hs @@ -0,0 +1,72 @@ +{----------------------------------------------------------------- + + (c) 2008-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. + +--------------------------------------------------------------------} + +-- | PrettyPrint provides tools for colored output to the terminal +module PrettyPrint ( Color(..) + , Intensity(..) + , color_string + , color_string_int + ) where + + +-- imports +import Prelude + + +-- define colors +data Color = Black | Red | Green | Yellow | Blue | Magenta + | Cyan | White | Reset + deriving(Enum) + + +-- define intensity +data Intensity = Normal | Bold + deriving(Eq) + + +-- convert a color into the corresponding color code string +get_color_code :: Color -> String +get_color_code = show . fromEnum + + +-- convert an intensity to the corresponding color code string +get_intensity_code :: Intensity -> String +get_intensity_code x + | x == Normal = "22" + | x == Bold = "1" + | otherwise = "1" + + +-- convenience wrapper around color string for bold colors +color_string :: Color -> String -> String +color_string = color_string_int Bold + + +-- convert a standard string into one graphically rendered +-- allows customization of color and intensity +color_string_int :: Intensity -> Color -> String -> String +color_string_int intensity col str = + "\ESC[" + ++ (get_intensity_code intensity) + ++ ";3" + ++ (get_color_code col) + ++ "m" + ++ str + ++ "\ESC[0;m" diff --git a/src/RpnCalc.hs b/src/RpnCalc.hs new file mode 100644 --- /dev/null +++ b/src/RpnCalc.hs @@ -0,0 +1,64 @@ +{----------------------------------------------------------------- + + (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. + +--------------------------------------------------------------------} + +-- | 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 + + +-- imports +import qualified Data.Map as M +import Prelude + +-- local imports +import GenericModel +import RpnData + + +-- | computes an expressions based on an rpn stack +-- 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 + + 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 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 item = item:ys diff --git a/src/RpnData.hs b/src/RpnData.hs new file mode 100644 --- /dev/null +++ b/src/RpnData.hs @@ -0,0 +1,52 @@ +{----------------------------------------------------------------- + + (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. + +--------------------------------------------------------------------} + +-- | RpnCalc defines the data structures and a calculator engine +-- for computing mathematical expressions that have been parsed +-- into reverse polish notations +module RpnData ( RpnItem(..) + , RpnStack + ) where + + +-- imports +import Prelude + + +-- | RpnItem describes all items that can be present in our +-- rpn stack +data RpnItem = Time + | Number Double + | Variable String + | UnaFunc (Double -> Double) + | BinFunc (Double -> Double -> Double) + +type RpnStack = [RpnItem] + + +instance Show RpnItem where + show (Time) = "TIME" + show (Number x) = show x + show (BinFunc x) = show (x 1 2) + show (UnaFunc x) = show (x 1) + show (Variable x) = x + + + diff --git a/src/RpnParser.hs b/src/RpnParser.hs new file mode 100644 --- /dev/null +++ b/src/RpnParser.hs @@ -0,0 +1,185 @@ +{----------------------------------------------------------------- + + (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. + +--------------------------------------------------------------------} + +-- | module parsing an infix math equation into an compute stack +-- in reverse polish notation (rpn) +module RpnParser ( parse_infix_to_rpn ) where + +-- imports +import Control.Monad +import Prelude + +-- local imports +import ExtraFunctions +import GenericModel +import RpnData +import TokenParser + +-- import Debug.Trace + +-- | 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) + "infix math expression" + + +-- | parser for expressions chained via "+" +add_term :: CharParser ModelState RpnStack +add_term = concat . insert_adds <$> + sub_term `sepBy` (reservedOp "+") + "addition term" + + where + insert_adds [] = [] + insert_adds (x:xs) = x:foldr (\y a -> y:[BinFunc (+)]:a) [] xs + + +-- | parser for expressions chained via "-" +sub_term :: CharParser ModelState RpnStack +sub_term = concat . insert_subs <$> + div_term `sepBy` (reservedOp "-") + "subtraction term" + + where + insert_subs [] = [] + insert_subs (x:xs) = x:foldr (\y a -> y:[BinFunc (-)]:a) [] xs + + +-- | parser for expressions chained via "*" +div_term :: CharParser ModelState RpnStack +div_term = concat . insert_divs <$> + mul_term `sepBy` (reservedOp "/") + "division term" + + where + insert_divs [] = [] + insert_divs (x:xs) = x:foldr (\y a -> y:[BinFunc (/)]:a) [] xs + + +-- | parser for expressions chained via "/" +mul_term :: CharParser ModelState RpnStack +mul_term = concat . insert_muls <$> + exp_term `sepBy` (reservedOp "*") + "product term" + + where + insert_muls [] = [] + 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 "^") + "exponent" + + where + insert_exps [] = [] + 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 = try parse_single_number -- need try due to the unary "-" + <|> try signed_parenthesis -- (otherwise we get stuck) + <|> parse_functions + <|> parse_variable + "token or variable" + + +-- | parse all operations of type (Double -> Double) +-- we currently know about +parse_functions :: CharParser ModelState RpnStack +parse_functions = (msum $ extract_ops builtinFunctions) + "builtin unary function" + where + extract_ops = foldr (\(x,y) acc -> + ((reserved x *> execute y):acc)) [] + execute op = ( (insert_func op) <$> + ( parens add_term + <|> parse_single_number + <|> parse_variable + )) + "function parsing" + + insert_func _ [] = [] + 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 = push_parens <$> parse_sign <*> parens add_term + "signed parenthesis" + + where + 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) + "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 +parse_number = converter <$> parse_sign <*> naturalOrFloat + "signed integer or double" + where + converter sign val = case val of + Left i -> sign * (fromInteger i) + 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) ) + "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) + "variable" + where + -- in case of a unary minus we also push the necessary + -- multiplication by (-1) onto the stack + push sign x = + if sign >= 0 + then + if x == "TIME" + then [Time] + else [Variable x] + else + if x == "TIME" + then [Time,Number (-1.0),BinFunc (*)] + else [Variable x, Number (-1.0), BinFunc(*)] + diff --git a/src/TokenParser.hs b/src/TokenParser.hs new file mode 100644 --- /dev/null +++ b/src/TokenParser.hs @@ -0,0 +1,184 @@ +{----------------------------------------------------------------- + + (c) 2008-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. + +--------------------------------------------------------------------} + +-- | functionality related to parsing tokens +module TokenParser ( module Control.Applicative + , module Text.ParserCombinators.Parsec + , braces + , builtinFunctions + , comma + , charLiteral + , float + , integer + , parens + , keywords + , lexer + , naturalOrFloat + , operators + , reservedOp + , reserved + , semi + , stringLiteral + , symbol + , whiteSpace + ) where + + +-- imports +import Control.Applicative +import Control.Monad (ap, MonadPlus (..)) +import Prelude +import Text.ParserCombinators.Parsec hiding (many,optional, (<|>)) +import qualified Text.ParserCombinators.Parsec.Token as PT +import Text.ParserCombinators.Parsec.Language (haskellDef + , opLetter + , reservedOpNames + , reservedNames ) + +-- local imports +import ExtraFunctions + + +{- Definitions for Applicative Parsec instance -} + +-- | Applicative instance for Monad +instance Applicative (GenParser s a) where + pure = return + (<*>) = ap + + +-- |Alternative instance for MonadPlus +instance Alternative (GenParser s a) where + empty = mzero + (<|>) = mplus + + + +{- set up the Token Parser -} + +-- | builtin functions of the form (Double -> Double) +builtinFunctions :: [(String, Double -> Double)] +builtinFunctions = [ ("sqrt",sqrt) + , ("exp",exp) + , ("log",log) + , ("log2", logBase 2) + , ("log10", logBase 10) + , ("sin",sin) + , ("cos",cos) + , ("tan",tan) + , ("asin", asin) + , ("acos", acos) + , ("atan", atan) + , ("sinh", sinh) + , ("cosh", cosh) + , ("tanh", tanh) + , ("asinh", sinh) + , ("acosh", cosh) + , ("atanh", atanh) + , ("erf", erf) + , ("erfc", erfc) + , ("abs", abs) + ] + + +-- | all other keywords that are not regular functions +keywords :: [String] +keywords = [ "def", "molecules", "reactions", "time", "outputIter" + , "nil", "outputFreq", "outputFile", "systemVol" + ] + +operators :: [String] +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 Integer +integer :: CharParser st Integer +integer = PT.integer lexer + + +-- | token parser for Char +stringLiteral :: CharParser st String +stringLiteral = PT.stringLiteral lexer + + +-- | token parser for Char +charLiteral :: CharParser st Char +charLiteral = PT.charLiteral lexer + + +-- | token parser for Double +float :: CharParser st Double +float = PT.float lexer + + +-- | token parser for Either Integer Double +naturalOrFloat :: CharParser st (Either Integer Double) +naturalOrFloat = PT.naturalOrFloat lexer + + +-- | token parser for keywords +reservedOp :: String -> CharParser st () +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 semicolon +semi :: CharParser st String +semi = PT.semi lexer + + +-- | token parser for comma +comma :: CharParser st String +comma = PT.comma lexer + + +-- | token parser for symbol +symbol :: String -> CharParser st String +symbol = PT.symbol lexer + diff --git a/src/simgi.hs b/src/simgi.hs new file mode 100644 --- /dev/null +++ b/src/simgi.hs @@ -0,0 +1,88 @@ +{----------------------------------------------------------------- + + (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. + +--------------------------------------------------------------------} + +-- | main gsim driver +module Main where + +-- imports +import Prelude +import System.IO +import System.Random +import System.Environment + +-- local imports +import CommandLine +import Engine +import InputCheck +import InputParser +import Messages +import TokenParser + +-- import Debug.Trace + + +-- | main +main :: IO () +main = + + -- process command line arguments + getArgs >>= process_commandline + >>= \((_,_),files) -> + + -- reject anything but a single input file + if length files /= 1 + then usage + else + + -- read input file and extract content + readFile (head files) + >>= \content -> + + -- parse input file + case runParser input_parser initialModelState "" content of + Left er -> putStrLn (show er) + Right parsedState -> + case check_input parsedState of + Left err -> putStrLn err + Right _ -> newStdGen + >>= \gen -> + + -- set up simuation + let + rands = randomRs (0,1) gen :: [Double] + initialOutput = create_initial_output parsedState + initialState = create_initial_state parsedState + rands initialOutput + totalTime = maxTime parsedState + dataDumpIter = maxIter parsedState + outFile = outfileName parsedState + in + + -- open output file + openFile outFile WriteMode + >>= \handle -> startup_message + + -- ready to run the simulation + >> gillespie_driver handle totalTime dataDumpIter + initialState + + -- close output file + >> hClose handle + diff --git a/test/Makefile b/test/Makefile new file mode 100644 --- /dev/null +++ b/test/Makefile @@ -0,0 +1,15 @@ + + +default: check + + +.PHONY: check clean +check: + make -C rpnstack_test + make -C reversible_test + + +clean: + make -C rpnstack_test clean + make -C reversible_test clean + diff --git a/test/reversible_test/Makefile b/test/reversible_test/Makefile new file mode 100644 --- /dev/null +++ b/test/reversible_test/Makefile @@ -0,0 +1,19 @@ + +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 diff --git a/test/reversible_test/average.hs b/test/reversible_test/average.hs new file mode 100644 --- /dev/null +++ b/test/reversible_test/average.hs @@ -0,0 +1,34 @@ +{----------------------------------------------------------------- + + (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) diff --git a/test/reversible_test/check_deviation.hs b/test/reversible_test/check_deviation.hs new file mode 100644 --- /dev/null +++ b/test/reversible_test/check_deviation.hs @@ -0,0 +1,45 @@ +{----------------------------------------------------------------- + + (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 + diff --git a/test/reversible_test/reversible.sgl b/test/reversible_test/reversible.sgl new file mode 100644 --- /dev/null +++ b/test/reversible_test/reversible.sgl @@ -0,0 +1,25 @@ +{----------------------------------------------------- + + 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 diff --git a/test/reversible_test/run_test.sh b/test/reversible_test/run_test.sh new file mode 100644 --- /dev/null +++ b/test/reversible_test/run_test.sh @@ -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} 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} diff --git a/test/rpnstack_test/Makefile b/test/rpnstack_test/Makefile new file mode 100644 --- /dev/null +++ b/test/rpnstack_test/Makefile @@ -0,0 +1,15 @@ + +SRC_LOCATION=../../src + + +default: check + + +.PHONY: check clean +check: + ghc -i$(SRC_LOCATION) --make RpnStackTest.hs + @./RpnStackTest + + +clean: + rm -f *.hi *.o RpnStackTest diff --git a/test/rpnstack_test/RpnStackTest.hs b/test/rpnstack_test/RpnStackTest.hs new file mode 100644 --- /dev/null +++ b/test/rpnstack_test/RpnStackTest.hs @@ -0,0 +1,227 @@ +{----------------------------------------------------------------- + + (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 "" "" "" +