diff --git a/.github/workflows/main-testing.yml b/.github/workflows/main-testing.yml new file mode 100644 index 00000000..4e076f71 --- /dev/null +++ b/.github/workflows/main-testing.yml @@ -0,0 +1,108 @@ +name: main-validation + +on: + push: + branches: [ master, issue37] + pull_request: + branches: [ master ] + +env: + # Customize the CMake build type here (Release, Debug, RelWithDebInfo, etc.) + BUILD_TYPE: Release + +jobs: + build: + # The CMake configure and build commands are platform agnostic and should work equally + # well on Windows or Mac. You can convert this to a matrix build if you need + # cross-platform coverage. + # See: https://docs.github.com/en/free-pro-team@latest/actions/learn-github-actions/managing-complex-workflows#using-a-build-matrix + runs-on: ${{ matrix.os }} + + strategy: + fail-fast: false + matrix: + os: [ubuntu-22.04, macos-latest, windows-latest] + + steps: + - uses: actions/checkout@v3 + - uses: seanmiddleditch/gha-setup-ninja@v3 + + # - name: Build on windows and run check + # if: matrix.os == 'windows-latest' + # shell: pwsh + # run: | + # mkdir build + # cd build + # cmake -G "Ninja" .. + # ninja + + # # Dynamically locate the MinGW bin directory + # $mingw_bin = (Get-Command gcc).Source | Split-Path + + # # Copy required DLLs into the build directory + # Copy-Item "$mingw_bin\libgcc_s_seh-1.dll" -Destination ${{ github.workspace }}\build -ErrorAction Stop + # Copy-Item "$mingw_bin\libstdc++-6.dll" -Destination ${{ github.workspace }}\build -ErrorAction Stop + # Copy-Item "$mingw_bin\libwinpthread-1.dll" -Destination ${{ github.workspace }}\build -ErrorAction Stop + + # # Run the executable to check it's working + # ./NFsim.exe -h + + + - name: Build on windows and run check (static linking) + if: matrix.os == 'windows-latest' + shell: pwsh + run: | + mkdir build + cd build + cmake -G "Ninja" -DCMAKE_EXE_LINKER_FLAGS="-static-libgcc -static-libstdc++ -static" .. + ninja + ./NFsim.exe -h + + + - name: Configure CMake + if: matrix.os != 'windows-latest' + # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using a single-configuration generator such as make. + # See https://cmake.org/cmake/help/latest/variable/CMAKE_BUILD_TYPE.html?highlight=cmake_build_type + run: cmake -B ${{github.workspace}}/build -DCMAKE_BUILD_TYPE=${{env.BUILD_TYPE}} + + - name: Build unix + if: matrix.os != 'windows-latest' + # Build your program with the given configuration + run: cmake --build ${{github.workspace}}/build --config ${{env.BUILD_TYPE}} + + - name: Set up Python + uses: actions/setup-python@v2 + with: + python-version: '3.11' + - name: Cache pip + uses: actions/cache@v3 + with: + path: ~/.cache/pip + key: ${{ runner.os }}-pip-${{ hashFiles('requirements.txt') }} + restore-keys: | + ${{ runner.os }}-pip- + ${{ runner.os }}- + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install -r validate/requirements.txt + - name: Validation + run: python validate/validate.py validate/. + - name: Prepare unix bundle + if: matrix.os != 'windows-latest' + run: | + cd + mkdir ${{github.workspace}}/bundle_${{ matrix.os }} + cp ${{github.workspace}}/build/NFsim ${{github.workspace}}/bundle_${{ matrix.os }}/. + - name: Prepare windows bundle + if: matrix.os == 'windows-latest' + run: | + cd + mkdir ${{github.workspace}}/bundle_${{ matrix.os }} + cp ${{github.workspace}}/build/NFsim.exe ${{github.workspace}}/bundle_${{ matrix.os }}/. + #cp ${{github.workspace}}/build/*.dll ${{github.workspace}}/bundle_${{ matrix.os }}/. + - name: Archive compiled NFsim + uses: actions/upload-artifact@v4 + with: + name: NFsim_bin_${{ matrix.os }} + path: ${{github.workspace}}/bundle_${{ matrix.os }} diff --git a/.gitignore b/.gitignore new file mode 100644 index 00000000..6848116a --- /dev/null +++ b/.gitignore @@ -0,0 +1,14 @@ +build +.vscode +build_test_files +validate/basicModels/*.*dat +validate/basicModels/*.net +validate/basicModels/*.xml +validate/basicModels/*.tsv +validate/basicModels/*.json +test/Issue37/test2.nfevent.json +test/Issue37/test.nfevent.json +test/Issue37/out.gdat +test/Issue37/issue37_test2.species +test/Issue37/issue37_test2.gdat +bin/NFsim diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 00000000..2a9d7f71 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,71 @@ +cmake_minimum_required(VERSION 2.8.9) + +set(CMAKE_OSX_ARCHITECTURES "x86_64;arm64" CACHE STRING "Build universal binary" FORCE) + +project(NFsim) + +#find_package(XMLRPC REQUIRED abyss-server c++2) + + +#if (NOT XMLRPC_FOUND) +# message (FATAL_ERROR "XMLRPC-C not found") +#endif (NOT XMLRPC_FOUND) + +set(CMAKE_BUILD_TYPE Release) + +set(SUB_DIRS + src/nauty24 + src/NFutil + src/NFutil/MTrand + src/NFtest/transcription + src/NFtest/tlbr + src/NFtest/simple_system + src/NFtest/agentcell/cell + src/NFtest/agentcell + src/NFscheduler + src/NFreactions/transformations + src/NFreactions/reactions + src/NFreactions/reactantLists + src/NFreactions/mappings + src/NFreactions + src/NFoutput + src/NFinput + src/NFinput/TinyXML + src/NFinput/json + src/NFfunction/muParser + src/NFfunction + src/NFcore + src/NFcore/reactionSelector + src/NFcore/moleculeLists +) +add_definitions("-Wno-deprecated-declarations") + + +IF(CMAKE_SYSTEM_NAME MATCHES "Windows") + SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11") +# SET(CMAKE_FIND_LIBRARY_SUFFIXES ".a") +# SET(BUILD_SHARED_LIBRARIES OFF) +# SET(CMAKE_EXE_LINKER_FLAGS "-static") +ELSE(CMAKE_SYSTEM_NAME MATCHES "Windows") + set_source_files_properties( src/NFreactions/reactions/DORreaction.cpp PROPERTIES COMPILE_FLAGS -std=c++11 ) + set_source_files_properties( src/NFcore/molecule.cpp PROPERTIES COMPILE_FLAGS -std=c++11 ) + set_source_files_properties( src/NFcore/templateMolecule.cpp PROPERTIES COMPILE_FLAGS -std=c++11 ) + set_source_files_properties( src/NFinput/NFinput.cpp PROPERTIES COMPILE_FLAGS -std=c++11 ) +ENDIF(CMAKE_SYSTEM_NAME MATCHES "Windows") + + +#include_directories(${CMAKE_CURRENT_SOURCE_DIR} ${XMLRPC_INCLUDE_DIRS} ${SUB_DIRS}) +include_directories(${CMAKE_CURRENT_SOURCE_DIR} include ${SUB_DIRS}) + + +file(GLOB_RECURSE src_files "${CMAKE_CURRENT_SOURCE_DIR}/src/*cpp") +file(GLOB_RECURSE c_src_files "${CMAKE_CURRENT_SOURCE_DIR}/src/*c") + + +set(SRC_FILES + ${src_files} + ${c_src_files} +) + + +add_executable(${PROJECT_NAME} ${SRC_FILES} ) diff --git a/CMakeLists.x86.txt b/CMakeLists.x86.txt new file mode 100644 index 00000000..b30e74a8 --- /dev/null +++ b/CMakeLists.x86.txt @@ -0,0 +1,58 @@ +cmake_minimum_required(VERSION 2.8.9) + +project(NFsim) + +#find_package(XMLRPC REQUIRED abyss-server c++2) + +SET(CMAKE_C_COMPILER C:/cygwin/bin/gcc.exe) +SET(CMAKE_CXX_COMPILER C:/cygwin/bin/g++.exe) + + +#if (NOT XMLRPC_FOUND) +# message (FATAL_ERROR "XMLRPC-C not found") +#endif (NOT XMLRPC_FOUND) + +set(CMAKE_BUILD_TYPE Release) + +set(SUB_DIRS + src/nauty24 + src/NFutil + src/NFutil/MTrand + src/NFtest/transcription + src/NFtest/tlbr + src/NFtest/simple_system + src/NFtest/agentcell/cell + src/NFtest/agentcell + src/NFscheduler + src/NFreactions/transformations + src/NFreactions/reactions + src/NFreactions/reactantLists + src/NFreactions/mappings + src/NFreactions + src/NFoutput + src/NFinput + src/NFinput/TinyXML + src/NFinput/json + src/NFfunction/muParser + src/NFfunction + src/NFcore + src/NFcore/reactionSelector + src/NFcore/moleculeLists +) +add_definitions("-Wno-deprecated-declarations") +set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11 -DTIXML_USE_STL ") + +include_directories(${CMAKE_CURRENT_SOURCE_DIR} include ${SUB_DIRS}) + + +file(GLOB_RECURSE src_files "${CMAKE_CURRENT_SOURCE_DIR}/src/*cpp") +file(GLOB_RECURSE c_src_files "${CMAKE_CURRENT_SOURCE_DIR}/src/*c") + + +set(SRC_FILES + ${src_files} + ${c_src_files} +) + + +add_executable(${PROJECT_NAME} ${SRC_FILES} ) diff --git a/LICENSE.txt b/LICENSE.txt index 818433ec..c74c937d 100644 --- a/LICENSE.txt +++ b/LICENSE.txt @@ -1,674 +1,21 @@ - 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 -. +MIT License + +Copyright (c) 2016 Michael W. Sneddon, James R. Faeder, Thierry Emonet + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/NFsim_manual_v1.11.pdf b/NFsim_manual_v1.11.pdf deleted file mode 100644 index 17725d80..00000000 Binary files a/NFsim_manual_v1.11.pdf and /dev/null differ diff --git a/NFsim_manual_v1.12.pdf b/NFsim_manual_v1.12.pdf new file mode 100644 index 00000000..c333cec0 Binary files /dev/null and b/NFsim_manual_v1.12.pdf differ diff --git a/README.md b/README.md new file mode 100644 index 00000000..bc1e0d33 --- /dev/null +++ b/README.md @@ -0,0 +1,267 @@ +# NFsim - the network free stochastic simulator + +[![NFsim build status](https://github.com/RuleWorld/nfsim/workflows/main-validation/badge.svg)](https://github.com/RuleWorld/nfsim/actions) + + Coverity Scan Build Status + + + +- michael w. sneddon +- justin s. hogg +- jose-juan tapia +- james r. faeder +- thierry emonet + +Yale University +University of Pittsburgh +funded by the National Science Foundation + +## Overview + +NFsim is a free, open-source, biochemical reaction simulator designed to handle +systems that have a large or even infinite number of possible molecular +interactions or states. NFsim also has advanced and flexible options for +simulating coarse-grained representations of complex nonlinear reaction +mechanisms. + +NFsim is ideal for modeling polymerization, aggregation, and cooperative +reactions that cannot be handled with traditional stochastic or ODE simulators. +Models are specified in the BioNetGen Langauge, providing a powerful model +building environment. + +If you just want to download and use NFsim, you should simply download a +preconfigured packaged release from http://emonet.biology.yale.edu/nfsim. If +you want to hack on the code or make contributions, please create a fork and +submit pull requests to the dev branch. + +If you use NFsim for your research or work, please cite NFsim as: Sneddon MW, +Faeder JR & Emonet T. Efficient modeling, simulation and coarse-graining of +biological complexity with NFsim. Nature Methods,(2011) 8(2):177-83. + +## Repository Contents + +NFsim is released under the MIT License. See LICENSE.txt for more details about +redistribution restrictions. + +For help with running NFsim, see the user manual, NFsim\_manual\_[version].pdf, +and open the example model "simple\_system.bngl". + +Source code is in the "src" directory. Example models are in the "models" +directory, with README files. BioNetGen and the ODE and SSA solvers used by +BioNetGen are in the BNG and Network2 directories. + +Enjoy your new network-free world! + +## Building NFsim + +### Linux and OSX + +Make sure you have a recent version of CMake installed and run the following at +a terminal: + + mkdir build + cd build + cmake .. + make + +### Windows + +Make sure you have a recent version of CMake, Cygwin (or MinGW), and Ninja +installed. Then run the following at a terminal: + + mkdir build + cd build + cmake -G "Ninja" + ninja + +If you've built it with Cygwin and you want to move or package up the +executable, you'll need to copy the the following DLLs along with it: + +* cygwin1.dll +* cygstdc++-6.dll +* cygz.dll +* cyggcc_s-seh-1.dll. + +## Release Notes + +### v1.14.3 April, 2025 + +Changed compilation flags for MacOS to produce universal binary the runs on both arm64 and intel macs. + +### v1.14.2 May, 2024 + +(a) Bugfix: Fixed a segmentation fault resulting from having a species start with zero concentration. Two simple models were added to the validation process to catch this error. +(b) Changes to the validation process were made. Previously, 15 trajectories would be generated for each model before testing the differences of NFsim and SSA versus NFsim and ODE. Now only one trajectory would be tested and more would generate upon faliure until a model failed 15 times in a row. + +### v1.14.1 May, 2023 + +(a) Bugfix: Added a missing "\[" when no operations are present + +### v1.14.0 February, 2023 + +Relabeling the version number to follow previous scheme to avoid further confusion. + +(a) Feature: New output format, `.nfevent.json`. JSON based format for all events that happen during a NFsim simulation. This format replaces the older reaction format and can be used with the `-rxnlog` command line argument. The use of `.nfevent.json` file extension is highly recommended since a JSON schema for the format is available. +(b) Bugfix: Fixed an issue where `-connect` option would break with models that have remove operations. + +### v1.2.2 April, 2022 + +(a) Bugfix: Disabled maxcputime option +(b) Bugfix: NFsim now checks for compartments and quits if it finds a compartments block +(c) Bugifx: Leftover debug statement was removed. + +### v1.2.1 Feb, 2022 + +(a) Bugfix release. An accidentally leftover debug statement was removed. + +### v1.2.0 Jan, 2022 + +(a) Added a new built-in function called TFUN that can pull values from a file given a counter observable in a model. +(b) A new option is added to infer connectivity between reaction rules. This option allows the user to infer "connected" rules before a model is ran and each time a rule fires, only connected rules are checked for updates, instead of every possible rule in the system. + +### v1.12.1 Aug, 2016 + +(a) Bugfix release. Addresses an error dealing with local function and species +labels. The error dealt with the way mappingSets where created and passed to +the localFunction evaluation. A test case (v19.bngl) was added that addresses +this case. + +### v1.12 Dec, 2015 + +(a) Changes to how molecule instances are mapped to BasicRxn's +(BasicRxnClass::tryToAdd()). It was possible for certain kinds of rules that +the mappingSets were not updated correctly because the head molecule matched a +reactant pattern before and after a reaction event BUT the mapping was +different after the firing. In such cases, the mappingSet was not updated +properly. The logic was changed to fix this problem. (b) Further changes to +how molecule instances are mapped to a ReactionClass object (BasicRxn and +DORReactions). In particular it is often the case that graph symmetry leads to +a complex being able to map to a ReactionClass multiple times. Symmetry +considerations were being made on the reaction center but not on the context +components which led to an undercounting of the number of times a pattern agent +could match a rule instance in some edge cases (see v17.bngl in the validation +suite). This led to incorrect results or even NFSim crashes. NOTE: the changes +made in points (a) and (b) may cause some models to execute less efficiently. +(c) Fixed index bound checking in MoleculeType::getComponentStateName(). (d) +Updated the validation models to reflect current BNGL formatting standards. +Added a few new validation models that address the bugfixes included in this +version. This update also includes a Python version of the validation script. + +### v1.11 Oct, 2012 + +(a) Molecules without components may be treated as population variables +rather than individual agents. This feature is useful for reducing memory +requirements when a simple molecule has very large population. A molecule +type may be flagged for treatment as a population variable by using the +"population" keyword following the molecule type definition in the BNGL +model file. See section 8.i in the documentation for further detail. +(b) Added a new Reaction Class called "FunctionProduct" that permits local +functions defined on two reactants. The local rate law must have the form +f(x)\*g(y), where x and y are tags on two distinct reactants. See section +7.c in the documentation for complete information. +(c) Fixed some problems evaluating complex-scoped local functions (CSLF). +CSLF were not updated properly after reactions that split complexes or +deleted molecules. As in v1.10, complex-scoped local functions are +enabled by default. A new command-line switch, -nocslf, has been added +which disables complex-scoped evaluation. +(d) Improved efficiency for matching patterns with connected-to syntax +when the connected-to component does not have reaction center. This may +be especially notable in models with large complexes, e.g. polymerization. + +### v1.10 Aug, 2011 + +(a) Command line parser now detects arguments that are not properly preceeded +by a dash, and generates a warning. (b) Includes a check when creating +template molecules that throws an error when users attempt to use Null or Trash +in reactant or observable patterns (anything that requires the creation of a +Template molecule). (c) fixed a bug introduced in v1.09 whereby a site was +allowed to bind to itself, for instance, in a dimerization rxn. (d) Support +for creating a new molecule bound to an existing molecule, as in a rule like +A(a) -> A(a!1).A(a!1). Existing code that implemented this feature did not +function properly with the check for null conditions before reactions were +fired. (e) Fixed bug in template molecule when clearing molecules after a +connected-to syntax search. In some cases, not all molecules were being +cleared, giving rise to situations where adding one observable created dangling +matches which affected the results of other observables. (f) NFsim is now +packaged with Network3, an updated version of the run\_network code to execute +ODE and SSA simulations. Network3 allows global functions in BNGL models among +other release features given here: +http://bionetgen.org/index.php/Release\_Notes. Note that Network3 does not +support On-The-Fly Stochastic Simulation (you will have to recompile Network2 +to use this feature). (g) Mac 32bit is no longer supported by NFsim, but you +can make executables for older Macs by recompiling the code on your own +machine. See the manual for instructions. + +### v1.09 Apr, 2011 + +(a) NFsim now allows the mixing of integers and strings as component labels, +although if numbers and strings are mixed, all labels are parsed as strings, +NOT integers. Therefore, PLUS and MINUS keywords cannot be used if mixing +integer states and string states, and a warning will be generated if a state is +set to PLUS. One can only use PLUS or MINUS when ALL states are an integer +value greater than zero. This new behavior was needed to handle BNGL files that +used the convention of ~P specifying phosphorylated, and ~0 specifying +unphosphorylated. (b) Fixed bug whereby if verbose option was turned on +without specifying an output file location, no output would be generated. Now, +output to a gdat file will be generated in these cases. (c) Use of 'ss' input +argument to 'saveSpecies', which prints a a list of all species at the end of a +simulation, is now handled. This feature was implemented to allow future +support in BNG by restarting an NFsim simulation after it ends, which can be +done by parsing the output species list together with a BNGL model file. This +feature still has to be tested in BNG, and will likely be fully documented in +v1.10. The 'ss' flag writes the file to either system_name_nf.species, or a +file designated by the user. (d) fixed memory leak in TemplateMolecules that +caused memory / performance issues with molecules having multiple identical +sites and a high degree of aggregation. (e) fixed csv error, where when the csv +flag is used, the header line is not comma delimited. (f) nfsim now supports +intra-molecular binding. Previously these events were rejected as null events. + +### v1.08 Dec, 2010 + +With the new TotalRate keyword, users are now able to specify whether or not to +use the microscopic (default) interpretation or macroscopic (TotalRate) +interpretation of rate laws. Now, NFsim convention matches BNG. Previously, +NFsim interpreted all rates as microscopic except for global functions, which +were interpreted as macroscopic. This is now also explained in the user manual. +Example models for the flagellar motor and oscillating gene expression have +been updated correspondingly so that they still produce the same results as in +the NFsim paper. + +Users also now have the option of outputting gdat files in a comma delimited +format (csv), which makes parsing the output file easier in some circumstances, +using the flag "-csv". Additionally, a bug in the parameter scanning script was +fixed that caused the script to crash when scanning a model that includes the +local function syntax. + +### v1.07 Nov, 2010 + +A series of updates to the code were made in this release. (1) RNF files that +are not found produce an error message. Previously, no error message was given +and execution proceeded as if the RNF flag was not given. (2) Input flags to +NFsim can be given in the original format with a single dash (as in ./NFsim +-logo), or with the more "linuxy" style double dash (as in ./NFsim --logo). +(3) The parameter scan script had problems when parsing BNGL files with local +functions due to the '%' character. This is now fixed. (4) The universal +traversal limit is automatically set to be the size of the largest pattern in +the system, which can be overridden by passing the -utl flag. This allows +users who are unfamiliar with this speedup to still take advantage of it to a +certain extent. (5) the -rtag flag was added that allows NFsim to produce +output whenever a particular reaction, given by the -rtag flag, is given. This +allows, for instance, users to track the fates of single particles exactly +without using the comprehensive molecule output feature. (6) The above changes +are documented in an updated user manual. + +### v1.06 Sept 28, 2010 + +Added scripts for running NFsim from Matlab, parameter scanning, and basic +parameter estimation. The manual is also updated to reflect these changes. +However, the precompiled executables of NFsim remain unchanged from v1.05, so +running them will give the old version number unless you recompile the code on +your own computer. Also, models that were used to compare the performance of +NFsim to DYNSTOC, RuleMonkey, and Kappa are now included with a readme file +under: models/performance\_test\_models. + +### v1.052 + +First publicly released stable build diff --git a/README.txt b/README.txt deleted file mode 100644 index 44ad11c1..00000000 --- a/README.txt +++ /dev/null @@ -1,188 +0,0 @@ - - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% % -% @@ @ @@@@@ % -% @ @ @ @ % -% @ @ @ @@@@ ___ % -% @ @ @ @ /__ | |\ /| % -% @ @@ @ ___\ | | v | % -% % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -NFsim - the network free stochastic simulator, v1.11 - -michael w. sneddon -justin s. hogg -james r. faeder -thierry emonet - -Yale University -University of Pittsburgh -funded by the National Science Foundation - - -################################################################################ - -NFsim is a free, open-source, biochemical reaction simulator designed to handle systems -that have a large or even infinite number of possible molecular interactions or states. -NFsim also has advanced and flexible options for simulating coarse-grained representations -of complex nonlinear reaction mechanisms. - -NFsim is ideal for modeling polymerization, aggregation, and cooperative reactions that -cannot be handled with traditional stochastic or ODE simulators. Models are specified in -the BioNetGen Langauge, providing a powerful model building environment. - -If you just want to download and use NFsim, you should simply download a preconfigured -packaged release from http://emonet.biology.yale.edu/nfsim. If you want to hack on the -code or make contributions, please create a fork and submit pull requests to the dev -branch. - -If you use NFsim for your research or work, please cite NFsim as: -Sneddon MW, Faeder JR & Emonet T. Efficient modeling, simulation and -coarse-graining of biological complexity with NFsim. Nature Methods,(2011) 8(2):177-83. - -################################################################################ - -NFsim is released under the GNU General Public License. See LICENSE.txt for -more details about redistribution restrictions. The git repository is now hosted -on github at: https://github.com/msneddon/nfsim - -For help with running NFsim, see the user manual, NFsim_manual_[version].pdf, -and open the example model "simple_system.bngl". - -Executable files for Windows, Mac and Linux are in the "bin" directory. Source -code and makefiles for NFsim are in the NFcode directory. Example models are -in the "models" directory, with README files. BioNetGen and the ODE and SSA -solvers used by BioNetGen are in the BNG and Network2 directories. Finally, -a suite of helpful analysis and other modeling tools can be found in the -"NFtools" directory. - -Enjoy your new network-free world! - - - -################################################################################ - -Release Notes - -v1.11 Oct, 2012 - (a) Molecules without components may be treated as population variables - rather than individual agents. This feature is useful for reducing memory - requirements when a simple molecule has very large population. A molecule - type may be flagged for treatment as a population variable by using the - "population" keyword following the molecule type definition in the BNGL - model file. See section 8.i in the documentation for further detail. - (b) Added a new Reaction Class called "FunctionProduct" that permits local - functions defined on two reactants. The local rate law must have the form - f(x)*g(y), where x and y are tags on two distinct reactants. See section - 7.c in the documentation for complete information. - (c) Fixed some problems evaluating complex-scoped local functions (CSLF). - CSLF were not updated properly after reactions that split complexes or - deleted molecules. As in v1.10, complex-scoped local functions are - enabled by default. A new command-line switch, -nocslf, has been added - which disables complex-scoped evaluation. - (d) Improved efficiency for matching patterns with connected-to syntax - when the connected-to component does not have reaction center. This may - be especially notable in models with large complexes, e.g. polymerization. - -v1.10 Aug, 2011 - (a) Command line parser now detects arguments that are not properly - preceeded by a dash, and generates a warning. (b) Includes a check when - creating template molecules that throws an error when users attempt to - use Null or Trash in reactant or observable patterns (anything that - requires the creation of a Template molecule). (c) fixed a bug introduced - in v1.09 whereby a site was allowed to bind to itself, for instance, in a - dimerization rxn. (d) Support for creating a new molecule bound to an - existing molecule, as in a rule like A(a) -> A(a!1).A(a!1). Existing code - that implemented this feature did not function properly with the check - for null conditions before reactions were fired. (e) Fixed bug in template - molecule when clearing molecules after a connected-to syntax search. In - some cases, not all molecules were being cleared, giving rise to situations - where adding one observable created dangling matches which affected the - results of other observables. (f) NFsim is now packaged with Network3, an - updated version of the run_network code to execute ODE and SSA simulations. - Network3 allows global functions in BNGL models among other release features - given here: http://bionetgen.org/index.php/Release_Notes. Note that Network3 - does not support On-The-Fly Stochastic Simulation (you will have to recompile - Network2 to use this feature). (g) Mac 32bit is no longer supported by NFsim, - but you can make executables for older Macs by recompiling the code on your - own machine. See the manual for instructions. - - -v1.09 Apr, 2011 - (a) NFsim now allows the mixing of integers and strings as component - labels, although if numbers and strings are mixed, all labels are parsed - as strings, NOT integers. Therefore, PLUS and MINUS keywords cannot - be used if mixing integer states and string states, and a warning will - be generated if a state is set to PLUS. One can only use PLUS or MINUS - when ALL states are an integer value greater than zero. This new behavior - was needed to handle BNGL files that used the convention of ~P specifying - phosphorylated, and ~0 specifying unphosphorylated. - (b) Fixed bug whereby if verbose option was turned on without specifying - an output file location, no output would be generated. Now, output to - a gdat file will be generated in these cases. - (c) Use of 'ss' input argument to 'saveSpecies', which prints a a list of - all species at the end of a simulation, is now handled. This feature was - implemented to allow future support in BNG by restarting an NFsim simulation - after it ends, which can be done by parsing the output species list together - with a BNGL model file. This feature still has to be tested in BNG, and will - likely be fully documented in v1.10. The 'ss' flag writes the file to either - [system_name]_nf.species, or a file designated by the user. - (d) fixed memory leak in TemplateMolecules that caused memory / performance - issues with molecules having multiple identical sites and a high degree - of aggregation. (e) fixed csv error, where when the csv flag is used, the - header line is not comma delimited. (f) nfsim now supports intra-molecular - binding. Previously these events were rejected as null events. - - -v1.08 Dec, 2010 - With the new TotalRate keyword, users are now able to - specify whether or not to use the microscopic (default) interpretation - or macroscopic (TotalRate) interpretation of rate laws. Now, NFsim - convention matches BNG. Previously, NFsim interpreted all rates as - microscopic except for global functions, which were interpreted as - macroscopic. This is now also explained in the user manual. Example - models for the flagellar motor and oscillating gene expression have - been updated correspondingly so that they still produce the same results - as in the NFsim paper. - - Users also now have the option of outputting gdat files in a comma - delimited format (csv), which makes parsing the output file easier - in some circumstances, using the flag "-csv". Additionally, a bug in - the parameter scanning script was fixed that caused the script to - crash when scanning a model that includes the local function syntax. - - -v1.07 Nov, 2010 - A series of updates to the code were made in this - release. (1) RNF files that are not found produce an error message. - Previously, no error message was given and execution proceeded as if - the RNF flag was not given. (2) Input flags to NFsim can be given - in the original format with a single dash (as in ./NFsim -logo), or - with the more "linuxy" style double dash (as in ./NFsim --logo). - (3) The parameter scan script had problems when parsing BNGL files - with local functions due to the '%' character. This is now fixed. - (4) The universal traversal limit is automatically set to be the size - of the largest pattern in the system, which can be overridden by - passing the -utl flag. This allows users who are unfamiliar with - this speedup to still take advantage of it to a certain extent. (5) - the -rtag flag was added that allows NFsim to produce output whenever - a particular reaction, given by the -rtag flag, is given. This allows, - for instance, users to track the fates of single particles exactly - without using the comprehensive molecule output feature. (6) The above - changes are documented in an updated user manual. - - -v1.06 Sept 28, 2010 - added scripts for running NFsim from Matlab, parameter - scanning, and basic parameter estimation. The manual is also updated - to reflect these changes. However, the precompiled executables of NFsim - remain unchanged from v1.05, so running them will give the old version - number unless you recompile the code on your own computer. Also, models - that were used to compare the performance of NFsim to DYNSTOC, RuleMonkey, - and Kappa are now included with a readme file under: - models/performance_test_models. - -v1.052 First publicly released stable build - - - - diff --git a/bin/src/NFcore/moleculeLists/subdir.mk b/bin/src/NFcore/moleculeLists/subdir.mk deleted file mode 100644 index 5671bb0f..00000000 --- a/bin/src/NFcore/moleculeLists/subdir.mk +++ /dev/null @@ -1,24 +0,0 @@ -################################################################################ -# Automatically-generated file. Do not edit! -################################################################################ - -# Add inputs and outputs from these tool invocations to the build variables -CPP_SRCS += \ -../src/NFcore/moleculeLists/moleculeList.cpp - -OBJS += \ -./src/NFcore/moleculeLists/moleculeList.o - -CPP_DEPS += \ -./src/NFcore/moleculeLists/moleculeList.d - - -# Each subdirectory must supply rules for building sources it contributes -src/NFcore/moleculeLists/%.o: ../src/NFcore/moleculeLists/%.cpp - @echo 'Building file: $<' - @echo 'Invoking: GCC C++ Compiler' - g++ -O3 -Wall -c -fmessage-length=0 -MMD -MP -MF"$(@:%.o=%.d)" -MT"$(@:%.o=%.d)" -o "$@" "$<" - @echo 'Finished building: $<' - @echo ' ' - - diff --git a/bin/src/NFcore/reactionSelector/subdir.mk b/bin/src/NFcore/reactionSelector/subdir.mk deleted file mode 100644 index 9d44dbbc..00000000 --- a/bin/src/NFcore/reactionSelector/subdir.mk +++ /dev/null @@ -1,27 +0,0 @@ -################################################################################ -# Automatically-generated file. Do not edit! -################################################################################ - -# Add inputs and outputs from these tool invocations to the build variables -CPP_SRCS += \ -../src/NFcore/reactionSelector/directSelector.cpp \ -../src/NFcore/reactionSelector/logClassSelector.cpp - -OBJS += \ -./src/NFcore/reactionSelector/directSelector.o \ -./src/NFcore/reactionSelector/logClassSelector.o - -CPP_DEPS += \ -./src/NFcore/reactionSelector/directSelector.d \ -./src/NFcore/reactionSelector/logClassSelector.d - - -# Each subdirectory must supply rules for building sources it contributes -src/NFcore/reactionSelector/%.o: ../src/NFcore/reactionSelector/%.cpp - @echo 'Building file: $<' - @echo 'Invoking: GCC C++ Compiler' - g++ -O3 -Wall -c -fmessage-length=0 -MMD -MP -MF"$(@:%.o=%.d)" -MT"$(@:%.o=%.d)" -o "$@" "$<" - @echo 'Finished building: $<' - @echo ' ' - - diff --git a/bin/src/NFcore/subdir.mk b/bin/src/NFcore/subdir.mk deleted file mode 100644 index 8710438b..00000000 --- a/bin/src/NFcore/subdir.mk +++ /dev/null @@ -1,45 +0,0 @@ -################################################################################ -# Automatically-generated file. Do not edit! -################################################################################ - -# Add inputs and outputs from these tool invocations to the build variables -CPP_SRCS += \ -../src/NFcore/complex.cpp \ -../src/NFcore/complexList.cpp \ -../src/NFcore/molecule.cpp \ -../src/NFcore/moleculeType.cpp \ -../src/NFcore/observable.cpp \ -../src/NFcore/reactionClass.cpp \ -../src/NFcore/system.cpp \ -../src/NFcore/templateMolecule.cpp - -OBJS += \ -./src/NFcore/complex.o \ -./src/NFcore/complexList.o \ -./src/NFcore/molecule.o \ -./src/NFcore/moleculeType.o \ -./src/NFcore/observable.o \ -./src/NFcore/reactionClass.o \ -./src/NFcore/system.o \ -./src/NFcore/templateMolecule.o - -CPP_DEPS += \ -./src/NFcore/complex.d \ -./src/NFcore/complexList.d \ -./src/NFcore/molecule.d \ -./src/NFcore/moleculeType.d \ -./src/NFcore/observable.d \ -./src/NFcore/reactionClass.d \ -./src/NFcore/system.d \ -./src/NFcore/templateMolecule.d - - -# Each subdirectory must supply rules for building sources it contributes -src/NFcore/%.o: ../src/NFcore/%.cpp - @echo 'Building file: $<' - @echo 'Invoking: GCC C++ Compiler' - g++ -O3 -Wall -c -fmessage-length=0 -MMD -MP -MF"$(@:%.o=%.d)" -MT"$(@:%.o=%.d)" -o "$@" "$<" - @echo 'Finished building: $<' - @echo ' ' - - diff --git a/bin/src/NFfunction/muParser/subdir.mk b/bin/src/NFfunction/muParser/subdir.mk deleted file mode 100644 index 48e31f46..00000000 --- a/bin/src/NFfunction/muParser/subdir.mk +++ /dev/null @@ -1,45 +0,0 @@ -################################################################################ -# Automatically-generated file. Do not edit! -################################################################################ - -# Add inputs and outputs from these tool invocations to the build variables -CPP_SRCS += \ -../src/NFfunction/muParser/muParser.cpp \ -../src/NFfunction/muParser/muParserBase.cpp \ -../src/NFfunction/muParser/muParserBytecode.cpp \ -../src/NFfunction/muParser/muParserCallback.cpp \ -../src/NFfunction/muParser/muParserComplex.cpp \ -../src/NFfunction/muParser/muParserError.cpp \ -../src/NFfunction/muParser/muParserInt.cpp \ -../src/NFfunction/muParser/muParserTokenReader.cpp - -OBJS += \ -./src/NFfunction/muParser/muParser.o \ -./src/NFfunction/muParser/muParserBase.o \ -./src/NFfunction/muParser/muParserBytecode.o \ -./src/NFfunction/muParser/muParserCallback.o \ -./src/NFfunction/muParser/muParserComplex.o \ -./src/NFfunction/muParser/muParserError.o \ -./src/NFfunction/muParser/muParserInt.o \ -./src/NFfunction/muParser/muParserTokenReader.o - -CPP_DEPS += \ -./src/NFfunction/muParser/muParser.d \ -./src/NFfunction/muParser/muParserBase.d \ -./src/NFfunction/muParser/muParserBytecode.d \ -./src/NFfunction/muParser/muParserCallback.d \ -./src/NFfunction/muParser/muParserComplex.d \ -./src/NFfunction/muParser/muParserError.d \ -./src/NFfunction/muParser/muParserInt.d \ -./src/NFfunction/muParser/muParserTokenReader.d - - -# Each subdirectory must supply rules for building sources it contributes -src/NFfunction/muParser/%.o: ../src/NFfunction/muParser/%.cpp - @echo 'Building file: $<' - @echo 'Invoking: GCC C++ Compiler' - g++ -O3 -Wall -c -fmessage-length=0 -MMD -MP -MF"$(@:%.o=%.d)" -MT"$(@:%.o=%.d)" -o "$@" "$<" - @echo 'Finished building: $<' - @echo ' ' - - diff --git a/bin/src/NFfunction/subdir.mk b/bin/src/NFfunction/subdir.mk deleted file mode 100644 index 8619ce05..00000000 --- a/bin/src/NFfunction/subdir.mk +++ /dev/null @@ -1,33 +0,0 @@ -################################################################################ -# Automatically-generated file. Do not edit! -################################################################################ - -# Add inputs and outputs from these tool invocations to the build variables -CPP_SRCS += \ -../src/NFfunction/compositeFunction.cpp \ -../src/NFfunction/funcParser.cpp \ -../src/NFfunction/function.cpp \ -../src/NFfunction/localFunction.cpp - -OBJS += \ -./src/NFfunction/compositeFunction.o \ -./src/NFfunction/funcParser.o \ -./src/NFfunction/function.o \ -./src/NFfunction/localFunction.o - -CPP_DEPS += \ -./src/NFfunction/compositeFunction.d \ -./src/NFfunction/funcParser.d \ -./src/NFfunction/function.d \ -./src/NFfunction/localFunction.d - - -# Each subdirectory must supply rules for building sources it contributes -src/NFfunction/%.o: ../src/NFfunction/%.cpp - @echo 'Building file: $<' - @echo 'Invoking: GCC C++ Compiler' - g++ -O3 -Wall -c -fmessage-length=0 -MMD -MP -MF"$(@:%.o=%.d)" -MT"$(@:%.o=%.d)" -o "$@" "$<" - @echo 'Finished building: $<' - @echo ' ' - - diff --git a/bin/src/NFinput/TinyXML/subdir.mk b/bin/src/NFinput/TinyXML/subdir.mk deleted file mode 100644 index c27d2618..00000000 --- a/bin/src/NFinput/TinyXML/subdir.mk +++ /dev/null @@ -1,33 +0,0 @@ -################################################################################ -# Automatically-generated file. Do not edit! -################################################################################ - -# Add inputs and outputs from these tool invocations to the build variables -CPP_SRCS += \ -../src/NFinput/TinyXML/tinystr.cpp \ -../src/NFinput/TinyXML/tinyxml.cpp \ -../src/NFinput/TinyXML/tinyxmlerror.cpp \ -../src/NFinput/TinyXML/tinyxmlparser.cpp - -OBJS += \ -./src/NFinput/TinyXML/tinystr.o \ -./src/NFinput/TinyXML/tinyxml.o \ -./src/NFinput/TinyXML/tinyxmlerror.o \ -./src/NFinput/TinyXML/tinyxmlparser.o - -CPP_DEPS += \ -./src/NFinput/TinyXML/tinystr.d \ -./src/NFinput/TinyXML/tinyxml.d \ -./src/NFinput/TinyXML/tinyxmlerror.d \ -./src/NFinput/TinyXML/tinyxmlparser.d - - -# Each subdirectory must supply rules for building sources it contributes -src/NFinput/TinyXML/%.o: ../src/NFinput/TinyXML/%.cpp - @echo 'Building file: $<' - @echo 'Invoking: GCC C++ Compiler' - g++ -O3 -Wall -c -fmessage-length=0 -MMD -MP -MF"$(@:%.o=%.d)" -MT"$(@:%.o=%.d)" -o "$@" "$<" - @echo 'Finished building: $<' - @echo ' ' - - diff --git a/bin/src/NFinput/subdir.mk b/bin/src/NFinput/subdir.mk deleted file mode 100644 index 36f5da54..00000000 --- a/bin/src/NFinput/subdir.mk +++ /dev/null @@ -1,39 +0,0 @@ -################################################################################ -# Automatically-generated file. Do not edit! -################################################################################ - -# Add inputs and outputs from these tool invocations to the build variables -CPP_SRCS += \ -../src/NFinput/NFinput.cpp \ -../src/NFinput/commandLineParser.cpp \ -../src/NFinput/parseFuncXML.cpp \ -../src/NFinput/parseSymRxns.cpp \ -../src/NFinput/rnfRunner.cpp \ -../src/NFinput/walk.cpp - -OBJS += \ -./src/NFinput/NFinput.o \ -./src/NFinput/commandLineParser.o \ -./src/NFinput/parseFuncXML.o \ -./src/NFinput/parseSymRxns.o \ -./src/NFinput/rnfRunner.o \ -./src/NFinput/walk.o - -CPP_DEPS += \ -./src/NFinput/NFinput.d \ -./src/NFinput/commandLineParser.d \ -./src/NFinput/parseFuncXML.d \ -./src/NFinput/parseSymRxns.d \ -./src/NFinput/rnfRunner.d \ -./src/NFinput/walk.d - - -# Each subdirectory must supply rules for building sources it contributes -src/NFinput/%.o: ../src/NFinput/%.cpp - @echo 'Building file: $<' - @echo 'Invoking: GCC C++ Compiler' - g++ -O3 -Wall -c -fmessage-length=0 -MMD -MP -MF"$(@:%.o=%.d)" -MT"$(@:%.o=%.d)" -o "$@" "$<" - @echo 'Finished building: $<' - @echo ' ' - - diff --git a/bin/src/NFoutput/subdir.mk b/bin/src/NFoutput/subdir.mk deleted file mode 100644 index e123f675..00000000 --- a/bin/src/NFoutput/subdir.mk +++ /dev/null @@ -1,24 +0,0 @@ -################################################################################ -# Automatically-generated file. Do not edit! -################################################################################ - -# Add inputs and outputs from these tool invocations to the build variables -CPP_SRCS += \ -../src/NFoutput/NFoutput.cpp - -OBJS += \ -./src/NFoutput/NFoutput.o - -CPP_DEPS += \ -./src/NFoutput/NFoutput.d - - -# Each subdirectory must supply rules for building sources it contributes -src/NFoutput/%.o: ../src/NFoutput/%.cpp - @echo 'Building file: $<' - @echo 'Invoking: GCC C++ Compiler' - g++ -O3 -Wall -c -fmessage-length=0 -MMD -MP -MF"$(@:%.o=%.d)" -MT"$(@:%.o=%.d)" -o "$@" "$<" - @echo 'Finished building: $<' - @echo ' ' - - diff --git a/bin/src/NFreactions/mappings/subdir.mk b/bin/src/NFreactions/mappings/subdir.mk deleted file mode 100644 index 320e5aed..00000000 --- a/bin/src/NFreactions/mappings/subdir.mk +++ /dev/null @@ -1,30 +0,0 @@ -################################################################################ -# Automatically-generated file. Do not edit! -################################################################################ - -# Add inputs and outputs from these tool invocations to the build variables -CPP_SRCS += \ -../src/NFreactions/mappings/mapping.cpp \ -../src/NFreactions/mappings/mappingGenerator.cpp \ -../src/NFreactions/mappings/mappingSet.cpp - -OBJS += \ -./src/NFreactions/mappings/mapping.o \ -./src/NFreactions/mappings/mappingGenerator.o \ -./src/NFreactions/mappings/mappingSet.o - -CPP_DEPS += \ -./src/NFreactions/mappings/mapping.d \ -./src/NFreactions/mappings/mappingGenerator.d \ -./src/NFreactions/mappings/mappingSet.d - - -# Each subdirectory must supply rules for building sources it contributes -src/NFreactions/mappings/%.o: ../src/NFreactions/mappings/%.cpp - @echo 'Building file: $<' - @echo 'Invoking: GCC C++ Compiler' - g++ -O3 -Wall -c -fmessage-length=0 -MMD -MP -MF"$(@:%.o=%.d)" -MT"$(@:%.o=%.d)" -o "$@" "$<" - @echo 'Finished building: $<' - @echo ' ' - - diff --git a/bin/src/NFreactions/reactantLists/subdir.mk b/bin/src/NFreactions/reactantLists/subdir.mk deleted file mode 100644 index 5ddf61ad..00000000 --- a/bin/src/NFreactions/reactantLists/subdir.mk +++ /dev/null @@ -1,27 +0,0 @@ -################################################################################ -# Automatically-generated file. Do not edit! -################################################################################ - -# Add inputs and outputs from these tool invocations to the build variables -CPP_SRCS += \ -../src/NFreactions/reactantLists/reactantList.cpp \ -../src/NFreactions/reactantLists/reactantTree.cpp - -OBJS += \ -./src/NFreactions/reactantLists/reactantList.o \ -./src/NFreactions/reactantLists/reactantTree.o - -CPP_DEPS += \ -./src/NFreactions/reactantLists/reactantList.d \ -./src/NFreactions/reactantLists/reactantTree.d - - -# Each subdirectory must supply rules for building sources it contributes -src/NFreactions/reactantLists/%.o: ../src/NFreactions/reactantLists/%.cpp - @echo 'Building file: $<' - @echo 'Invoking: GCC C++ Compiler' - g++ -O3 -Wall -c -fmessage-length=0 -MMD -MP -MF"$(@:%.o=%.d)" -MT"$(@:%.o=%.d)" -o "$@" "$<" - @echo 'Finished building: $<' - @echo ' ' - - diff --git a/bin/src/NFreactions/reactions/subdir.mk b/bin/src/NFreactions/reactions/subdir.mk deleted file mode 100644 index 9dca5e36..00000000 --- a/bin/src/NFreactions/reactions/subdir.mk +++ /dev/null @@ -1,27 +0,0 @@ -################################################################################ -# Automatically-generated file. Do not edit! -################################################################################ - -# Add inputs and outputs from these tool invocations to the build variables -CPP_SRCS += \ -../src/NFreactions/reactions/DORreaction.cpp \ -../src/NFreactions/reactions/reaction.cpp - -OBJS += \ -./src/NFreactions/reactions/DORreaction.o \ -./src/NFreactions/reactions/reaction.o - -CPP_DEPS += \ -./src/NFreactions/reactions/DORreaction.d \ -./src/NFreactions/reactions/reaction.d - - -# Each subdirectory must supply rules for building sources it contributes -src/NFreactions/reactions/%.o: ../src/NFreactions/reactions/%.cpp - @echo 'Building file: $<' - @echo 'Invoking: GCC C++ Compiler' - g++ -O3 -Wall -c -fmessage-length=0 -MMD -MP -MF"$(@:%.o=%.d)" -MT"$(@:%.o=%.d)" -o "$@" "$<" - @echo 'Finished building: $<' - @echo ' ' - - diff --git a/bin/src/NFreactions/subdir.mk b/bin/src/NFreactions/subdir.mk deleted file mode 100644 index 1a8258a5..00000000 --- a/bin/src/NFreactions/subdir.mk +++ /dev/null @@ -1,24 +0,0 @@ -################################################################################ -# Automatically-generated file. Do not edit! -################################################################################ - -# Add inputs and outputs from these tool invocations to the build variables -CPP_SRCS += \ -../src/NFreactions/NFreactions.cpp - -OBJS += \ -./src/NFreactions/NFreactions.o - -CPP_DEPS += \ -./src/NFreactions/NFreactions.d - - -# Each subdirectory must supply rules for building sources it contributes -src/NFreactions/%.o: ../src/NFreactions/%.cpp - @echo 'Building file: $<' - @echo 'Invoking: GCC C++ Compiler' - g++ -O3 -Wall -c -fmessage-length=0 -MMD -MP -MF"$(@:%.o=%.d)" -MT"$(@:%.o=%.d)" -o "$@" "$<" - @echo 'Finished building: $<' - @echo ' ' - - diff --git a/bin/src/NFreactions/transformations/subdir.mk b/bin/src/NFreactions/transformations/subdir.mk deleted file mode 100644 index c4ea0a57..00000000 --- a/bin/src/NFreactions/transformations/subdir.mk +++ /dev/null @@ -1,33 +0,0 @@ -################################################################################ -# Automatically-generated file. Do not edit! -################################################################################ - -# Add inputs and outputs from these tool invocations to the build variables -CPP_SRCS += \ -../src/NFreactions/transformations/moleculeCreator.cpp \ -../src/NFreactions/transformations/speciesCreator.cpp \ -../src/NFreactions/transformations/transformation.cpp \ -../src/NFreactions/transformations/transformationSet.cpp - -OBJS += \ -./src/NFreactions/transformations/moleculeCreator.o \ -./src/NFreactions/transformations/speciesCreator.o \ -./src/NFreactions/transformations/transformation.o \ -./src/NFreactions/transformations/transformationSet.o - -CPP_DEPS += \ -./src/NFreactions/transformations/moleculeCreator.d \ -./src/NFreactions/transformations/speciesCreator.d \ -./src/NFreactions/transformations/transformation.d \ -./src/NFreactions/transformations/transformationSet.d - - -# Each subdirectory must supply rules for building sources it contributes -src/NFreactions/transformations/%.o: ../src/NFreactions/transformations/%.cpp - @echo 'Building file: $<' - @echo 'Invoking: GCC C++ Compiler' - g++ -O3 -Wall -c -fmessage-length=0 -MMD -MP -MF"$(@:%.o=%.d)" -MT"$(@:%.o=%.d)" -o "$@" "$<" - @echo 'Finished building: $<' - @echo ' ' - - diff --git a/bin/src/NFscheduler/subdir.mk b/bin/src/NFscheduler/subdir.mk deleted file mode 100644 index 5c6b15ac..00000000 --- a/bin/src/NFscheduler/subdir.mk +++ /dev/null @@ -1,27 +0,0 @@ -################################################################################ -# Automatically-generated file. Do not edit! -################################################################################ - -# Add inputs and outputs from these tool invocations to the build variables -CPP_SRCS += \ -../src/NFscheduler/NFstream.cpp \ -../src/NFscheduler/Scheduler.cpp - -OBJS += \ -./src/NFscheduler/NFstream.o \ -./src/NFscheduler/Scheduler.o - -CPP_DEPS += \ -./src/NFscheduler/NFstream.d \ -./src/NFscheduler/Scheduler.d - - -# Each subdirectory must supply rules for building sources it contributes -src/NFscheduler/%.o: ../src/NFscheduler/%.cpp - @echo 'Building file: $<' - @echo 'Invoking: GCC C++ Compiler' - g++ -O3 -Wall -c -fmessage-length=0 -MMD -MP -MF"$(@:%.o=%.d)" -MT"$(@:%.o=%.d)" -o "$@" "$<" - @echo 'Finished building: $<' - @echo ' ' - - diff --git a/bin/src/NFtest/agentcell/cell/subdir.mk b/bin/src/NFtest/agentcell/cell/subdir.mk deleted file mode 100644 index df30d208..00000000 --- a/bin/src/NFtest/agentcell/cell/subdir.mk +++ /dev/null @@ -1,30 +0,0 @@ -################################################################################ -# Automatically-generated file. Do not edit! -################################################################################ - -# Add inputs and outputs from these tool invocations to the build variables -CPP_SRCS += \ -../src/NFtest/agentcell/cell/cell.cpp \ -../src/NFtest/agentcell/cell/environment.cpp \ -../src/NFtest/agentcell/cell/util.cpp - -OBJS += \ -./src/NFtest/agentcell/cell/cell.o \ -./src/NFtest/agentcell/cell/environment.o \ -./src/NFtest/agentcell/cell/util.o - -CPP_DEPS += \ -./src/NFtest/agentcell/cell/cell.d \ -./src/NFtest/agentcell/cell/environment.d \ -./src/NFtest/agentcell/cell/util.d - - -# Each subdirectory must supply rules for building sources it contributes -src/NFtest/agentcell/cell/%.o: ../src/NFtest/agentcell/cell/%.cpp - @echo 'Building file: $<' - @echo 'Invoking: GCC C++ Compiler' - g++ -O3 -Wall -c -fmessage-length=0 -MMD -MP -MF"$(@:%.o=%.d)" -MT"$(@:%.o=%.d)" -o "$@" "$<" - @echo 'Finished building: $<' - @echo ' ' - - diff --git a/bin/src/NFtest/agentcell/subdir.mk b/bin/src/NFtest/agentcell/subdir.mk deleted file mode 100644 index b77ec6ff..00000000 --- a/bin/src/NFtest/agentcell/subdir.mk +++ /dev/null @@ -1,24 +0,0 @@ -################################################################################ -# Automatically-generated file. Do not edit! -################################################################################ - -# Add inputs and outputs from these tool invocations to the build variables -CPP_SRCS += \ -../src/NFtest/agentcell/agentcell.cpp - -OBJS += \ -./src/NFtest/agentcell/agentcell.o - -CPP_DEPS += \ -./src/NFtest/agentcell/agentcell.d - - -# Each subdirectory must supply rules for building sources it contributes -src/NFtest/agentcell/%.o: ../src/NFtest/agentcell/%.cpp - @echo 'Building file: $<' - @echo 'Invoking: GCC C++ Compiler' - g++ -O3 -Wall -c -fmessage-length=0 -MMD -MP -MF"$(@:%.o=%.d)" -MT"$(@:%.o=%.d)" -o "$@" "$<" - @echo 'Finished building: $<' - @echo ' ' - - diff --git a/bin/src/NFtest/simple_system/subdir.mk b/bin/src/NFtest/simple_system/subdir.mk deleted file mode 100644 index ef5ebd30..00000000 --- a/bin/src/NFtest/simple_system/subdir.mk +++ /dev/null @@ -1,24 +0,0 @@ -################################################################################ -# Automatically-generated file. Do not edit! -################################################################################ - -# Add inputs and outputs from these tool invocations to the build variables -CPP_SRCS += \ -../src/NFtest/simple_system/simple_system.cpp - -OBJS += \ -./src/NFtest/simple_system/simple_system.o - -CPP_DEPS += \ -./src/NFtest/simple_system/simple_system.d - - -# Each subdirectory must supply rules for building sources it contributes -src/NFtest/simple_system/%.o: ../src/NFtest/simple_system/%.cpp - @echo 'Building file: $<' - @echo 'Invoking: GCC C++ Compiler' - g++ -O3 -Wall -c -fmessage-length=0 -MMD -MP -MF"$(@:%.o=%.d)" -MT"$(@:%.o=%.d)" -o "$@" "$<" - @echo 'Finished building: $<' - @echo ' ' - - diff --git a/bin/src/NFtest/tlbr/subdir.mk b/bin/src/NFtest/tlbr/subdir.mk deleted file mode 100644 index 966a27f6..00000000 --- a/bin/src/NFtest/tlbr/subdir.mk +++ /dev/null @@ -1,24 +0,0 @@ -################################################################################ -# Automatically-generated file. Do not edit! -################################################################################ - -# Add inputs and outputs from these tool invocations to the build variables -CPP_SRCS += \ -../src/NFtest/tlbr/tlbr.cpp - -OBJS += \ -./src/NFtest/tlbr/tlbr.o - -CPP_DEPS += \ -./src/NFtest/tlbr/tlbr.d - - -# Each subdirectory must supply rules for building sources it contributes -src/NFtest/tlbr/%.o: ../src/NFtest/tlbr/%.cpp - @echo 'Building file: $<' - @echo 'Invoking: GCC C++ Compiler' - g++ -O3 -Wall -c -fmessage-length=0 -MMD -MP -MF"$(@:%.o=%.d)" -MT"$(@:%.o=%.d)" -o "$@" "$<" - @echo 'Finished building: $<' - @echo ' ' - - diff --git a/bin/src/NFtest/transcription/subdir.mk b/bin/src/NFtest/transcription/subdir.mk deleted file mode 100644 index 1b0aa6c6..00000000 --- a/bin/src/NFtest/transcription/subdir.mk +++ /dev/null @@ -1,24 +0,0 @@ -################################################################################ -# Automatically-generated file. Do not edit! -################################################################################ - -# Add inputs and outputs from these tool invocations to the build variables -CPP_SRCS += \ -../src/NFtest/transcription/transcription.cpp - -OBJS += \ -./src/NFtest/transcription/transcription.o - -CPP_DEPS += \ -./src/NFtest/transcription/transcription.d - - -# Each subdirectory must supply rules for building sources it contributes -src/NFtest/transcription/%.o: ../src/NFtest/transcription/%.cpp - @echo 'Building file: $<' - @echo 'Invoking: GCC C++ Compiler' - g++ -O3 -Wall -c -fmessage-length=0 -MMD -MP -MF"$(@:%.o=%.d)" -MT"$(@:%.o=%.d)" -o "$@" "$<" - @echo 'Finished building: $<' - @echo ' ' - - diff --git a/bin/src/NFutil/MTrand/subdir.mk b/bin/src/NFutil/MTrand/subdir.mk deleted file mode 100644 index b64c7726..00000000 --- a/bin/src/NFutil/MTrand/subdir.mk +++ /dev/null @@ -1,27 +0,0 @@ -################################################################################ -# Automatically-generated file. Do not edit! -################################################################################ - -# Add inputs and outputs from these tool invocations to the build variables -CPP_SRCS += \ -../src/NFutil/MTrand/mtrand.cpp \ -../src/NFutil/MTrand/mttest.cpp - -OBJS += \ -./src/NFutil/MTrand/mtrand.o \ -./src/NFutil/MTrand/mttest.o - -CPP_DEPS += \ -./src/NFutil/MTrand/mtrand.d \ -./src/NFutil/MTrand/mttest.d - - -# Each subdirectory must supply rules for building sources it contributes -src/NFutil/MTrand/%.o: ../src/NFutil/MTrand/%.cpp - @echo 'Building file: $<' - @echo 'Invoking: GCC C++ Compiler' - g++ -O3 -Wall -c -fmessage-length=0 -MMD -MP -MF"$(@:%.o=%.d)" -MT"$(@:%.o=%.d)" -o "$@" "$<" - @echo 'Finished building: $<' - @echo ' ' - - diff --git a/bin/src/NFutil/subdir.mk b/bin/src/NFutil/subdir.mk deleted file mode 100644 index 10d82d65..00000000 --- a/bin/src/NFutil/subdir.mk +++ /dev/null @@ -1,30 +0,0 @@ -################################################################################ -# Automatically-generated file. Do not edit! -################################################################################ - -# Add inputs and outputs from these tool invocations to the build variables -CPP_SRCS += \ -../src/NFutil/conversion.cpp \ -../src/NFutil/random.cpp \ -../src/NFutil/stringOperations.cpp - -OBJS += \ -./src/NFutil/conversion.o \ -./src/NFutil/random.o \ -./src/NFutil/stringOperations.o - -CPP_DEPS += \ -./src/NFutil/conversion.d \ -./src/NFutil/random.d \ -./src/NFutil/stringOperations.d - - -# Each subdirectory must supply rules for building sources it contributes -src/NFutil/%.o: ../src/NFutil/%.cpp - @echo 'Building file: $<' - @echo 'Invoking: GCC C++ Compiler' - g++ -O3 -Wall -c -fmessage-length=0 -MMD -MP -MF"$(@:%.o=%.d)" -MT"$(@:%.o=%.d)" -o "$@" "$<" - @echo 'Finished building: $<' - @echo ' ' - - diff --git a/bin/src/nauty24/subdir.mk b/bin/src/nauty24/subdir.mk deleted file mode 100644 index eb6bd868..00000000 --- a/bin/src/nauty24/subdir.mk +++ /dev/null @@ -1,30 +0,0 @@ -################################################################################ -# Automatically-generated file. Do not edit! -################################################################################ - -# Add inputs and outputs from these tool invocations to the build variables -C_SRCS += \ -../src/nauty24/nausparse.c \ -../src/nauty24/nautil.c \ -../src/nauty24/nauty.c - -OBJS += \ -./src/nauty24/nausparse.o \ -./src/nauty24/nautil.o \ -./src/nauty24/nauty.o - -C_DEPS += \ -./src/nauty24/nausparse.d \ -./src/nauty24/nautil.d \ -./src/nauty24/nauty.d - - -# Each subdirectory must supply rules for building sources it contributes -src/nauty24/%.o: ../src/nauty24/%.c - @echo 'Building file: $<' - @echo 'Invoking: GCC C Compiler' - gcc -O3 -Wall -c -fmessage-length=0 -MMD -MP -MF"$(@:%.o=%.d)" -MT"$(@:%.o=%.d)" -o "$@" "$<" - @echo 'Finished building: $<' - @echo ' ' - - diff --git a/bin/src/subdir.mk b/bin/src/subdir.mk deleted file mode 100644 index 5857f89e..00000000 --- a/bin/src/subdir.mk +++ /dev/null @@ -1,24 +0,0 @@ -################################################################################ -# Automatically-generated file. Do not edit! -################################################################################ - -# Add inputs and outputs from these tool invocations to the build variables -CPP_SRCS += \ -../src/NFsim.cpp - -OBJS += \ -./src/NFsim.o - -CPP_DEPS += \ -./src/NFsim.d - - -# Each subdirectory must supply rules for building sources it contributes -src/%.o: ../src/%.cpp - @echo 'Building file: $<' - @echo 'Invoking: GCC C++ Compiler' - g++ -O3 -Wall -c -fmessage-length=0 -MMD -MP -MF"$(@:%.o=%.d)" -MT"$(@:%.o=%.d)" -o "$@" "$<" - @echo 'Finished building: $<' - @echo ' ' - - diff --git a/dist/changeFilename.py b/dist/changeFilename.py new file mode 100644 index 00000000..6b7da8ed --- /dev/null +++ b/dist/changeFilename.py @@ -0,0 +1,12 @@ +import shutil +import platform +import os + +version = '1.14.3' +destdir = os.path.join(version, '{0}-{1}'.format(platform.system(), platform.architecture()[0])) +os.makedirs(destdir) +print('moving NFSIm to {0}\n'.format(destdir)) +if platform.system() != 'Windows': + shutil.move('NFsim', os.path.join(destdir, 'NFsim')) +else: + shutil.move('NFsim.exe', os.path.join(destdir, 'NFsim.exe')) diff --git a/doc/NFsim_manual_v1.11.docx b/doc/NFsim_manual_v1.11.docx deleted file mode 100755 index 65a66c57..00000000 Binary files a/doc/NFsim_manual_v1.11.docx and /dev/null differ diff --git a/doc/NFsim_manual_v1.12.docx b/doc/NFsim_manual_v1.12.docx new file mode 100755 index 00000000..a3901742 Binary files /dev/null and b/doc/NFsim_manual_v1.12.docx differ diff --git a/make_bundles b/make_bundles new file mode 100644 index 00000000..36664654 --- /dev/null +++ b/make_bundles @@ -0,0 +1,42 @@ +#!/bin/bash + +# Set the version variable +VERSION="1_14_3" # Replace with the desired version + +# Loop through all zip files in the current directory +for zipfile in *.zip; do + # Extract the base name of the zip file (without extension) + BASENAME=$(basename "$zipfile" .zip) + + # Determine the OS based on the zip file name + if [[ "$BASENAME" == *win* ]]; then + OS="win" + elif [[ "$BASENAME" == *mac* ]]; then + OS="mac" + elif [[ "$BASENAME" == *ubuntu* ]]; then + OS="linux" + else + echo "Unknown OS for file: $zipfile" + continue + fi + + # Create a temporary directory for extraction + TEMP_DIR=$(mktemp -d) + + # Unzip the file into the temporary directory + unzip -q "$zipfile" -d "$TEMP_DIR" + + # Find the executable file and set its permissions to 755 + find "$TEMP_DIR" -type f -executable -exec chmod 755 {} \; + + # Create the new tar.gz archive + OUTPUT_FILE="NFsim_${VERSION}_${OS}.tgz" + tar -czf "$OUTPUT_FILE" -C "$TEMP_DIR" . + + # Clean up the temporary directory + rm -rf "$TEMP_DIR" + + echo "Processed $zipfile -> $OUTPUT_FILE" +done + +echo "All files processed." \ No newline at end of file diff --git a/nfevent.schema.json b/nfevent.schema.json new file mode 100644 index 00000000..81ebb691 --- /dev/null +++ b/nfevent.schema.json @@ -0,0 +1,152 @@ + +{ + "$schema": "https://json-schema.org/draft/2020-12/schema", + "$id": "https://raw.githubusercontent.com/RuleWorld/nfsim/ASinanSaglam/issue35/nf.schema.json", + "title": "NFsim simulation", + "description": "Output from a NFsim simulation", + "type": "object", + "properties": { + "simulation": { + "type": "object", + "properties": { + "info": { + "description": "Information regarding the simulation", + "type": "object", + "properties": { + "name": { + "description": "Name of our simulation", + "type": "string" + }, + "global_molecule_limit": { + "description": "Number of global molecules allowed for the simulation", + "type": "number" + }, + "number_of_molecule_types": { + "description": "Number of total molecule types defined", + "type": "number" + }, + "number_of_molecules": { + "description": "Number of total starting molecules", + "type": "number" + } + }, + "required": [ "name" ] + }, + "molecule_types": { + "description": "Molecule types in a simulation", + "type": "array", + "items": { "$ref": "#/$defs/molecule_type" }, + "minItems": 1, + "uniqueItems": true + }, + "initialState": { + "description": "Initial state of the simulation", + "type": "object", + "properties": { + "molTypes": { + "description": "A run-length encoded array of molecule type IDs for the starting state", + "type": "array", + "items": { "$ref": "#/$defs/molType" }, + "minItems": 0, + "uniqueItems": false + }, + "ops": { + "description": "A set of operations that need to be performed on the initial molecules to get the inital state", + "type": "array", + "items": { "$ref": "#/$defs/operation" }, + "minItems": 0, + "uniqueItems": false + } + }, + "required": [ "molTypes", "ops" ] + }, + "firings": { + "description": "All firings in the simulation", + "type": "array", + "items": { "$ref": "#/$defs/firing" }, + "minItems": 0, + "uniqueItems": false + } + }, + "required": [ "info", "molecule_types", "initialState", "firings" ] + } + }, + "$defs": { + "molecule_type": { + "description": "Description of each molecule type, their names, IDs and components", + "type": "object", + "properties": { + "name": { + "description": "Name of the molecule type", + "type": "string" + }, + "typeID": { + "description": "Internal ID used for the molecule type, every operation refers to this ID", + "type": "number" + }, + "components": { + "description": "Names of each component for the molecule type", + "type": "array", + "items": { "type": "string" }, + "minItems": 0 + }, + "componentStates": { + "description": "Names of each component states, by default the first state is used upon initialization", + "type": "array", + "minItems": 0, + "items": { + "type": "array", + "items": { "type": "string" }, + "minItems": 0 + } + } + }, + "required": [ "name", "typeID", "components", "componentStates" ] + }, + "molecule_array": { + "description": "This is an item in a run-length encoded array, first item is the molecule type index, second item is how many times it repeats", + "type": "array", + "prefixItems": [ + {"type": "number"}, + {"type": "number"} + ], + "items": false, + "uniqueItems": false + }, + "operation": { + "description": "Array fully describing an operation and what it operates on. Each operation has different inputs.", + "type": "array", + "prefixItems": [ + {"enum": ["Add", "Delete", "StateChange", "AddBond", "DeleteBond", "IncrementState", "DecrementState", "DecrementPopulation"] }, + {"type": "number"} + ], + "items": { "type": "number" }, + "minItems": 2, + "maxItems": 5, + "uniqueItems": false + }, + "firing": { + "description": "Full description of a single event in the simulation.", + "type": "object", + "properties": { + "props": { + "description": "The name of the reaction that fired, the global event number followed by the internal time it happened", + "type": "array", + "prefixItems": [ + {"type": "string"}, + {"type": "number"}, + {"type": "number"} + ], + "items": false, + "uniqueItems": false + }, + "ops": { + "description": "The set of operations that occured during the event.", + "type": "array", + "items": { "$ref": "#/$defs/operation" } + } + }, + "required": [ "props", "ops" ] + } + } +} \ No newline at end of file diff --git a/src/NFcore/NFcore.hh b/src/NFcore/NFcore.hh index b36316e9..c9f92dbd 100644 --- a/src/NFcore/NFcore.hh +++ b/src/NFcore/NFcore.hh @@ -54,6 +54,7 @@ namespace NFcore //Forward declarations to deal with cyclic dependencies class MapGenerator; class MappingSet; + class Mapping; class ReactantList; class TransformationSet; class MoleculeList; @@ -100,6 +101,39 @@ namespace NFcore class ReactionSelector; + //exception for the handling of local functions and mapping sets + class LocalFunctionException: public exception + { + public: + void setType1_Mol(vector * type1_Mol){ + this->type1_Mol = type1_Mol; + } + + vector* getType1_Mol() const{ + return type1_Mol; + } + + void setIndex(const int index){ + this->index = index; + } + + int getIndex() const{ + return index; + } + + virtual const char* what() const throw() + { + return "Species scope parameter exception"; + } + + private: + vector* type1_Mol; + int index; + + + + + }; //! Container to organize all system complexes. /*! @@ -211,11 +245,16 @@ namespace NFcore double getAverageGroupValue(string groupName, int valIndex); ReactionClass *getReaction(int rIndex) { return allReactions.at(rIndex); }; + vector getAllReactions () { return allReactions; }; + ReactionClass * getReactionByName(string name); MoleculeType * getMoleculeType(int mtIndex) { return allMoleculeTypes.at(mtIndex); }; MoleculeType * getMoleculeTypeByName(string name); int getNumOfMoleculeTypes() { return allMoleculeTypes.size(); }; Molecule * getMoleculeByUid(int uid); + // AS2023 - this is added to allow for the suppression of warnings if the + // molecule with the UID doesn't exist + Molecule * getMoleculeByUid(int uid, bool warn); int getNumOfMolecules(); //Functions used when setting up the system. Note that most of these methods @@ -239,6 +278,10 @@ namespace NFcore bool saveSpecies() { return saveSpecies(string(name+"_nf.species")); }; bool saveSpecies(string filename); + // AS2023 - this gets set up by rxnlog argument and enables the logging + // of each firing in the system + bool getReactionTrackingStatus() { return reactionTrackingEnabled; }; + void setReactionTrackingStatus(bool status) { reactionTrackingEnabled = status; }; LocalFunction * getLocalFunctionByName(string fName); //bool addFunctionReference(FunctionReference *fr); @@ -256,6 +299,16 @@ namespace NFcore /* tell the system where to ouptut results*/ void setOutputToBinary(); void registerOutputFileLocation(string filename); + /* reaction firings are output to this file + * if any reaction has tag flag set to 1 */ + void registerReactionFileLocation(string filename); + /* Connected reactions upon each reaction firing are written to this location */ + void registerConnectedRxnFileLocation(string filename); + /* Connected reactions for each reaction as calculated are written to this location */ + void registerListOfConnectedRxnFileLocation(string filename); + /* list of molecule types and reaction firing counts are stored in these files */ + void registerMoleculeTypeFileLocation(string filename); + void registerRxnListFileLocation(string filename); void setDumpOutputter(DumpSystem *ds); @@ -265,7 +318,7 @@ namespace NFcore void turnOffGlobalFuncOut() { this->outputGlobalFunctionValues=false; }; /* once all elements are added, we need to prepare and initialize for simulations */ - void tagReaction(int rID); + void tagReaction(unsigned int rID); @@ -278,6 +331,8 @@ namespace NFcore void outputAllObservableCounts(); void outputAllObservableCounts(double cSampleTime); void outputAllObservableCounts(double cSampleTime,int eventCounter); + void outputAllMoleculeTypes(); + void outputAllRxnFiringCounts(); int getNumOfSpeciesObs() const; Observable * getSpeciesObs(int index) const; @@ -346,6 +401,7 @@ namespace NFcore void turnOff_OnTheFlyObs(); void turnOnOutputEventCounter() { outputEventCounter=true; }; + int getGlobalEventCounter() { return globalEventCounter; }; void addParameter(string name,double value); double getParameter(string name); @@ -354,6 +410,9 @@ namespace NFcore void printAllParameters(); NFstream& getOutputFileStream(); + NFstream& getReactionFileStream(); + NFstream& getConnectedRxnFileStream(); + NFstream& getConnectedRxnListFileStream(); // NETGEN -- method to access allComplexes ComplexList & getAllComplexes( ) { return allComplexes; }; @@ -368,7 +427,74 @@ namespace NFcore file is generated. */ void turnOnCSVformat() { this->csvFormat = true; }; + /* + * turn tagged reaction flag on for output if any reactions are tagged + * Arvind Rasi Subramaniam + */ + void turnOnTagRxnOutput() { this->anyRxnTagged = true; }; + bool getAnyRxnTagged() { return anyRxnTagged; }; + + /* + * track connected reactions after each reaction firing for + * debugging. + * Arvind Rasi Subramaniam Nov 21, 2018 + */ + void setTrackConnected(bool value) { this->trackConnected = value; }; + bool getTrackConnected() { return this->trackConnected; }; + /* + * print connected reactions at the beginning of the simulation for + * checking that we represented our intuition correctly in the model. + * Arvind Rasi Subramaniam Nov 12, 2019 + */ + void setPrintConnected(bool value) { this->printConnected = value; }; + bool getPrintConnected() { return this->printConnected; }; + /* + * Write reaction and molecule numbers instead of names to reduce file size + * This flag will also output two tsv files of list of molecules and reactions. + * Arvind Rasi Subramaniam Feb 18, 2019 + */ + void setRxnNumberTrack(bool value) { this->trackRxnNumber = value; }; + bool getRxnNumberTrack() { return this->trackRxnNumber; }; + + /* Turn on inference and use of connectivity + * Arvind Rasi Subramaniam + */ + void useConnectivityFlag(bool connectivityFlag) {this->connectivityFlag = connectivityFlag;}; + bool getConnectivityFlag() {return connectivityFlag;}; + + void setMaxCpuTime(double time) { max_cpu_time = time; }; + + clock_t start,finish; + double current_cpu_time = 0; + + void setConnectedReactions(int rxn1, int rxn2) {connectedReactions[rxn1][rxn2] = true;}; + bool areReactionsConnected(int rxn1, int rxn2) {return connectedReactions[rxn1][rxn2];}; + + /* Track the last reaction firing time + * This is done to make sure that the reaction runs to completion + * @author: Arvind R. Subramaniam + * @date: 13 Nov 2019 + */ + void setLastRxnTime(double rxnTime) {this->lastRxnTime = rxnTime;}; + double getLastRxnTime() {return this->lastRxnTime;}; + + /* + * print new data + */ + void setOutputMoleculeTypes(bool setVal) { this->outputMoleculeTypesFile = setVal; }; + void setOutputRxnFiringCounts(bool setVal) { this->outputRxnFiringCountsFile = setVal; }; + + // AS2023 - buffer size setter/getter, in units of firings. This is + // here to allow the user to juggle the added burden of reaction logging + // between memory vs CPU. + int getLogBufferSize() { return this->log_buffer_size; }; + void setLogBufferSize(int bsize) { this->log_buffer_size = bsize; }; + + // AS2023 - Species log helper functions + void setSpeciesLog(string logstr) { this->speciesLog = logstr; }; + string getSpeciesLog() { return this->speciesLog; }; + protected: /////////////////////////////////////////////////////////////////////////// @@ -384,13 +510,28 @@ namespace NFcore bool outputGlobalFunctionValues; /*< set to true to output the value of all global functions at each output step */ int globalMoleculeLimit; /*< total number of any particular molecule that can be created, default=100,000 */ bool outputEventCounter; /*< set to true to output the cumulative number of events at each output step */ + bool anyRxnTagged; /*< sets whether any reaction is tagged for output when it fires */ + bool connectivityFlag; /* Whether to infer and use reaction connectivity for updating molecule rxn membership*/ + bool trackConnected; /* Whether to track connected reactions after each reaction firing. Useful for debugging */ + bool printConnected; /* Whether to print connected reactions at the beginning of the simulation. Useful for debugging */ + bool outputMoleculeTypesFile; /* Output molecule types (default: false) */ + bool outputRxnFiringCountsFile; /* Output reaction firing counts (default: false) */ + bool trackRxnNumber; /* Whether to track reaction numbers instead of names for minimizing file size */ + double lastRxnTime; /* Time when the last reaction was fired */ + bool reactionTrackingEnabled = false; /* tells if reaction tracking is on, for rxnlog argument */ int globalEventCounter; + string speciesLog; /* AS2023 - log string for initial species */ + /////////////////////////////////////////////////////////////////////////// // The container objects that maintain the core system configuration vector allMoleculeTypes; /*!< container of all MoleculeTypes in the simulation */ vector allReactions; /*!< container of all Reactions in the simulation */ + /* created by rasi to see if the reaction connectivity can be inferred at the beginning + * to avoid the time-intensive search during each simulation run. + */ + vector allTemplateMolecules; /*!< container of all Reactions in the simulation */ // NETGEN //vector allComplexes; /*!< container of all complexes in the simulation */ //queue nextAvailableComplex; /*!< queue tells us which complexes can be used next */ @@ -419,18 +560,26 @@ namespace NFcore double a_tot; /*< the sum of all a's (propensities) of all reactions */ double current_time; /*< keeps track of the simulation time */ ReactionClass * nextReaction; /*< keeps track of the next reaction to fire */ + // max CPU time for simulation + double max_cpu_time; /////////////////////////////////////////////////////////////////////////// // protected functions needed only by the system while running a simulation double get_A_tot() const { return a_tot; }; double recompute_A_tot(); double getNextRxn(); + double getMaxCpuTime() const { return max_cpu_time; }; /////////////////////////////////////////////////////////////////////////// // Neccessary variables and methods for outputting //ofstream outputFileStream; /* the stream to a file to write out the results */ NFstream outputFileStream; /* NFstream is a smart stream that uses ofstream or stringstream depending on whether NF_MPI is defined */ + NFstream reactionOutputFileStream; /* NFstream is a smart stream that uses ofstream or stringstream depending on whether NF_MPI is defined */ + NFstream connectedRxnFileStream; /* NFstream is a smart stream that uses ofstream or stringstream depending on whether NF_MPI is defined */ + NFstream connectedRxnListFileStream; /* NFstream is a smart stream that uses ofstream or stringstream depending on whether NF_MPI is defined */ + NFstream moleculeTypeFileStream; + NFstream rxnListFileStream; void outputGroupDataHeader(); @@ -442,8 +591,7 @@ namespace NFcore /////////////////////////////////////////////////////////////////////////// //random data structures and variables used for optimization.... - int **rxnIndexMap; /*!< maps reaction index values to a reaction, used for MoleculeTypes to - quickly lookup a reaction */ + int **rxnIndexMap; map paramMap; @@ -453,6 +601,11 @@ namespace NFcore //Data structure that performs the selection of the next reaction class ReactionSelector * selector; + // To look up connected reactions quickly + vector > connectedReactions; + + // AS2023 - sets the default log buffer size to 10000 firings. + int log_buffer_size = 10000; private: list molList; @@ -533,6 +686,8 @@ namespace NFcore string getComponentName(int cIndex) const; void getPossibleComponentStates(int cIndex, list &nameList); int getDefaultComponentState(int cIndex) const { return defaultCompState[cIndex]; }; + // AS2023 - we need the comp states in a vector + vector < vector < string > > getPossibleCompStates() {return possibleCompStates;}; int getCompIndexFromName(string cName) const; string getComponentStateName(int cIndex, int cValue); @@ -628,6 +783,11 @@ namespace NFcore /* updates a molecules membership (assumes molecule is of type this) */ void updateRxnMembership(Molecule * m); + /* Updates only molecule membership in connected reactions. + * The connected reactions are inferred at the simulation start. + * Arvind Rasi Subramaniam + */ + void updateConnectedRxnMembership(Molecule * m, ReactionClass * r); /* auto populate with default molecules */ void populateWithDefaultMolecules(int moleculeCount); @@ -780,6 +940,17 @@ namespace NFcore bool incrementPopulation(); bool decrementPopulation(); + // Not sure why these functions were not there, but I made the visitation variables protected + // I left the remaining public for now + // Arvind Rasi Subramaniam + bool getVisitedMolecule() const { return hasVisitedMolecule; } + void setVisitedMolecule(bool visit) { hasVisitedMolecule = visit; } + bool * hasVisitedBond; + TemplateMolecule *isMatchedTo; + + /* used when reevaluating local functions */ + bool hasEvaluatedMolecule; + /////////////////////////////////////////////////////////////////////// int getComponentState(int cIndex) const { return component[cIndex]; }; int getComponentIndexOfBond(int cIndex) const { return indexOfBond[cIndex]; }; @@ -839,14 +1010,20 @@ namespace NFcore /* static functions which bind and unbind two molecules */ static void bind(Molecule *m1, int cIndex1, Molecule *m2, int cIndex2); static void bind(Molecule *m1, string compName1, Molecule *m2, string compName2); - static void unbind(Molecule *m1, int bSiteIndex); - static void unbind(Molecule *m1, char * bSiteName); + // AS2023 - unbind now returns the indices of the molecule that's selected + // for the unbinding for recording/logging purposesz + static vector unbind(Molecule *m1, int bSiteIndex); + static vector unbind(Molecule *m1, char * bSiteName); /* functions needed to traverse a complex and get all components * which is important when we want to update reactions and complexes */ void traverseBondedNeighborhood(list &members, int traversalLimit); + // AS2023 - additional call sig to use with reaction firing logging + void traverseBondedNeighborhood(list &members, int traversalLimit, string &logstr); static void breadthFirstSearch(list &members, Molecule *m, int depth); + // AS2023 - additional call sig to use with reaction firing logging + static void breadthFirstSearch(list &members, Molecule *m, int depth, string &logstr); void depthFirstSearch(list &members); /* when we are ready to begin simulations, moleculeType calls this function @@ -855,7 +1032,7 @@ namespace NFcore /* function that tells this molecule that it changed states or bonds * and it should update its reaction membership */ - void updateRxnMembership(); + void updateRxnMembership(ReactionClass * r, bool useConnectivity); void removeFromObservables(); void addToObservables(); //void updateDORs(); @@ -882,6 +1059,11 @@ namespace NFcore /* print functions for debugging */ void printDetails(); void printDetails(ostream &o); + + /* print bond details for debugging and tracking state changes */ + void printBondDetails(); + void printBondDetails(NFstream &o); + void printBondDetails(ostream &o); static void printMoleculeList(list &members); static int getUniqueIdCount() { return uniqueIdCount; }; @@ -892,13 +1074,8 @@ namespace NFcore void setIsObs(int oIndex, int isObs) { isObservable[oIndex]=isObs; }; - /* used for traversing a molecule complex */ - bool hasVisitedMolecule; - bool * hasVisitedBond; - TemplateMolecule *isMatchedTo; - - /* used when reevaluating local functions */ - bool hasEvaluatedMolecule; + // /* used when reevaluating local functions */ + // bool hasEvaluatedMolecule; static const int NOSTATE = -1; @@ -906,21 +1083,10 @@ namespace NFcore static const int NOINDEX = -1; - //void addDependentUpdateMolecule(Molecule *m); - //void removeDependentUpdateMolecule(Molecule *m); - //void traverseBondedNeighborhoodForUpdate(list &members, int traversalLimit); - - - //Unnecessary functions now that were once used to - //mark a molecule that is not in the simulation - //bool isMolAlive() const { return !isDead; }; - //bool isMolDead() const { return isDead; }; - //void kill() { isDead = true; }; - //void create() { isDead = false; }; - protected: - + /* used for traversing a molecule complex */ + bool hasVisitedMolecule; bool isPrepared; bool isAliveInSim; @@ -931,7 +1097,6 @@ namespace NFcore int ID_unique; int listId; - static int uniqueIdCount; /* The type of this molecule */ @@ -978,7 +1143,6 @@ namespace NFcore static queue q; static queue d; static list ::iterator molIter; - //static list ::iterator molIter2; }; @@ -1031,12 +1195,19 @@ namespace NFcore ReactionClass(string name, double rate, string baseRateParameterName, TransformationSet *transformationSet, System *s); virtual ~ReactionClass(); + //! To get all reactant and product templates for inferring connectivity between reactions + //! @author Arvind Rasi Subramaniam + void setAllReactantAndProductTemplates(map reactants, + map products); + int getNumOfReactants() const { return n_reactants; }; string getName() const { return name; }; + int getFireCounter() const { return fireCounter; }; double getBaseRate() const { return baseRate; }; int getRxnType() const { return reactionType; }; + MoleculeType *getMoleculeTypeOfReactantTemplate(int pos) const; void setBaseRate(double newBaseRate,string newBaseRateName); void resetBaseRateFromSystemParamter(); @@ -1045,6 +1216,9 @@ namespace NFcore double get_a() const { return a; }; virtual void printDetails() const; void fire(double random_A_number); + // AS2023 - additional call sig to use with reaction firing tracking. The call + // will now return a log the event + string fire(double random_A_number, bool track); //For DOR reactions virtual void notifyRateFactorChange(Molecule * m, int reactantIndex, int rxnListIndex) = 0; @@ -1081,6 +1255,11 @@ namespace NFcore void setTotalRateFlag(bool totalRate) { totalRateFlag = totalRate; }; + /* Whether to use reaction connecitivity for updating molecule + * membership + * Arvind Rasi Subramaniam + */ + void setConnectivityFlag(bool flag) { useConnectivity = flag; }; // _NETGEN_ @@ -1089,6 +1268,21 @@ namespace NFcore bool tagged; + // Use by Arvind Rasi Subramaniam to speed up simulations + // by inferring connectivity beforehand + void appendConnectedRxn(ReactionClass * rxn); + void appendPreConnectedRxn(ReactionClass * rxn); + bool isReactionConnected(ReactionClass * rxn); + int getNumConnectedRxns() {return connectedReactions.size();}; + ReactionClass * getconnectedRxn(int rxn2_id) {return connectedReactions.at(rxn2_id);}; + + // Methods to identify connected reactions within NFsim + // Gateway method + void identifyConnectedReactions(); + + // Called from within Transformation Set to check connectivity + bool areMoleculeTypeAndComponentPresent(MoleculeType * mt, int cIndex); + bool isTemplateCompatible(TemplateMolecule * t); protected: virtual void pickMappingSets(double randNumber) const=0; @@ -1111,6 +1305,17 @@ namespace NFcore unsigned int traversalLimit; + /* Used for scanning all reactants and products to infer connectivity + * Arvind Rasi Subramaniam + */ + vector allReactantTemplates; + vector allProductTemplates; + /* Maintain a list of connected reactions whose reactant numbers + * might change upon firing this reaction. + * Arvind Rasi Subramaniam + */ + vector connectedReactions; + TemplateMolecule **reactantTemplates; TransformationSet * transformationSet; MappingSet **mappingSet; @@ -1150,6 +1355,11 @@ namespace NFcore */ vector symmetricMappingSet; bool comparisonResult; + /* whether to use reaction connectivity for updating molecule + * membership + * Arvind Rasi Subramaniam + */ + bool useConnectivity; }; diff --git a/src/NFcore/complex.cpp b/src/NFcore/complex.cpp index d27a72dd..9b7f0620 100644 --- a/src/NFcore/complex.cpp +++ b/src/NFcore/complex.cpp @@ -110,7 +110,7 @@ void Complex::mergeWithList(Complex * c) // move molecules in c to this complex c->refactorToNewComplex(this->ID_complex); - this->complexMembers.splice(complexMembers.end(),c->complexMembers); + this->complexMembers.splice(complexMembers.end(),c->complexMembers); (system->getAllComplexes()).notifyThatComplexIsAvailable(c->getComplexID()); } @@ -143,7 +143,7 @@ void Complex::updateComplexMembership(Molecule * m) unsetCanonical(); //Get list of things this molecule is still connected to - list members; + list members; m->traverseBondedNeighborhood(members, ReactionClass::NO_LIMIT); //counter++; @@ -167,19 +167,17 @@ void Complex::updateComplexMembership(Molecule * m) //cout<<" forming new complex: next available: " <getComplexID()<::iterator molIter; + list ::iterator molIter; for( molIter = members.begin(); molIter != members.end(); molIter++ ) { (*molIter)->moveToNewComplex(newComplex->getComplexID()); } //put our new complex elements into that complex - newComplex->complexMembers.splice(newComplex->complexMembers.end(),members); + newComplex->complexMembers.splice(newComplex->complexMembers.end(),members); //cout<<"size of list now: " << members.size() <ID_complex)); - - + complexMembers.remove_if(IsInWrongComplex(this->ID_complex)); //update new complex in reactions? diff --git a/src/NFcore/complexList.cpp b/src/NFcore/complexList.cpp index 3a13e1ac..c482b9ba 100644 --- a/src/NFcore/complexList.cpp +++ b/src/NFcore/complexList.cpp @@ -180,17 +180,15 @@ void ComplexList::outputMoleculeTypeCountPerComplex(MoleculeType *m) } - -// TODO: figure out how friend functions work!! // friend functions -template -NFstream& operator<<(NFstream& nfstream, const T& value) -{ - if (nfstream.useFile_) - nfstream.file_ << value; - else - nfstream.str_ << value; - - return nfstream; -} +// template +// NFstream& operator<<(NFstream& nfstream, const T& value) +// { +// if (nfstream.useFile_) +// nfstream.file_ << value; +// else +// nfstream.str_ << value; + +// return nfstream; +// } diff --git a/src/NFcore/molecule.cpp b/src/NFcore/molecule.cpp index e31ba972..36fe2ead 100644 --- a/src/NFcore/molecule.cpp +++ b/src/NFcore/molecule.cpp @@ -158,9 +158,14 @@ LocalFunction * Molecule::getLocalFunction(int localFunctionIndex) { -void Molecule::updateRxnMembership() +void Molecule::updateRxnMembership(ReactionClass * r, bool useConnectivity) { - parentMoleculeType->updateRxnMembership(this); + if (useConnectivity) { + parentMoleculeType->updateConnectedRxnMembership(this, r); + } + else { + parentMoleculeType->updateRxnMembership(this); + } } void Molecule::updateTypeIIFunctions() @@ -344,11 +349,14 @@ void Molecule::printDetails(ostream &o) o<<" components: "; for(int c=0; cgetComponentName(c) <<"="; o<getComponentStateName(c,component[c]); o<<"\tbond="; - if(bond[c]==NOBOND) o<<"empty"; + if(bond[c]==nullptr) o<<"empty"; else { o<getMoleculeTypeName()<<"_"<getUniqueID(); o<<"("<getMoleculeType()->getComponentName(this->indexOfBond[c])<<")"; @@ -368,12 +376,82 @@ void Molecule::printDetails(ostream &o) } } +/** + * Prints to screen, see printDetails(ostram ) + * @author Arvind Rasi Subramaniam + */ +void Molecule::printBondDetails() { + this->printBondDetails(cout); +} + +/** + * Print all bonded states and their details to output stream + * + * Iterates through all components of a molecule type and prints + * only those components with a bond including the binding partner + * and the site of attachment. + * + * Same function written separately for printing to cout and NFStream + * @param o - Stream to write to + * @author Arvind Rasi Subramaniam + */ +void Molecule::printBondDetails(ostream &o) { + int degree = 0; + o<< parentMoleculeType->getName() << "\t"<getComponentName(c) <<"="; +// o<getComponentStateName(c,component[c]); +// o<<",bond="; +// o<getMoleculeTypeName()<<"_"<getUniqueID(); +// o<<"("<getMoleculeType()->getComponentName(this->indexOfBond[c])<<")"; +// } +// } + o.flush(); +} + +/** + * Same as printBondDetails(ostream ) + * but writes to NFStream file instead of ostream + * @author Arvind Rasi Subramaniam + */ +void Molecule::printBondDetails(NFstream &o) +{ + if (parentMoleculeType->getSystem()->getRxnNumberTrack()) { + o << parentMoleculeType->getTypeID() << "\t" << ID_unique; + } else { + o << parentMoleculeType->getName() << "\t" << ID_unique; + } + if (parentMoleculeType->getSystem()->getTrackConnected()) { + o<<"\t"; + for(int c=0; cgetComponentName(c); + if (parentMoleculeType->getComponentStateName(c,component[c]) != "NO_STATE") { + o<< "-" << parentMoleculeType->getComponentStateName(c,component[c]); + } + o<<":"; + o<getMoleculeType()->getComponentName(this->indexOfBond[c]); + o<<"-"<getMoleculeTypeName()<<"_"<getUniqueID(); + } + } + } + o.flush(); +} + //Get the number of molecules this molecule is bonded to int Molecule::getDegree() { int degree = 0; for(int c=0; cbond[cIndex1]!=NOBOND || m2->bond[cIndex2]!=NOBOND) { + if(m1->bond[cIndex1]!=nullptr || m2->bond[cIndex2]!=nullptr) { cerr< Molecule::unbind(Molecule *m1, int cIndex) { //get the other molecule bound to this site //cout<<"I am here. "<getComplex()->printDetails(); + + // AS2023 - this now returns what is unbound in a vector + vector tpl; + tpl.push_back(m2->getUniqueID()); + tpl.push_back(cIndex2); + return tpl; } -void Molecule::unbind(Molecule *m1, char * compName) +vector Molecule::unbind(Molecule *m1, char * compName) { int cIndex = m1->getMoleculeType()->getCompIndexFromName(compName); - Molecule::unbind(m1,cIndex); + return Molecule::unbind(m1,cIndex); } @@ -531,7 +616,6 @@ void Molecule::unbind(Molecule *m1, char * compName) queue Molecule::q; queue Molecule::d; list ::iterator Molecule::molIter; - void Molecule::breadthFirstSearch(list &members, Molecule *m, int depth) { if(m==0) { @@ -595,11 +679,78 @@ void Molecule::breadthFirstSearch(list &members, Molecule *m, int d (*molIter)->hasVisitedMolecule=false; } +// AS2023 - alternative call sig for logging that includes a log string +void Molecule::breadthFirstSearch(list &members, Molecule *m, int depth, string &logstr) +{ + if(m==0) { + cerr<<"Error in Molecule::breadthFirstSearch, m is null.\n"; + cerr<<"Likely an internal error where a MappingSet is on a list and\n"; + cerr<<"is not actually mapped to any molecule!"; + exit(3); + } + + //Create the queues (for effeciency, now queues are a static attribute of Molecule...) + //queue q; + //queue d; + int currentDepth = 0; + + //cout<<"traversing on:"<printDetails(); + + //First add this molecule + q.push(m); + members.push_back(m); + d.push(currentDepth+1); + m->hasVisitedMolecule=true; + + //Look at children until the queue is empty + while(!q.empty()) + { + //Get the next parent to look at (currentMolecule) + Molecule *cM = q.front(); + currentDepth = d.front(); + q.pop(); + d.pop(); + + if (!logstr.empty()) { + logstr += " [\"Delete\"," + to_string(cM->getUniqueID()) + "],\n"; + } + + //Make sure the depth does not exceed the limit we want to search + if((depth!=ReactionClass::NO_LIMIT) && (currentDepth>=depth)) continue; + + //Loop through the bonds + int cMax = cM->numOfComponents; + for(int c=0; cgetComp + if(cM->isBindingSiteBonded(c)) + { + Molecule *neighbor = cM->getBondedMolecule(c); + //cout<<"looking at neighbor: "<printDetails(); + if(!neighbor->hasVisitedMolecule) + { + neighbor->hasVisitedMolecule=true; + members.push_back(neighbor); + q.push(neighbor); + d.push(currentDepth+1); + //cout<<"adding... to traversal list."<hasVisitedMolecule=false; +} + -// void Molecule::traverseBondedNeighborhood(list &members, int traversalLimit) { //always call breadth first search, it is a bit faster @@ -609,6 +760,15 @@ void Molecule::traverseBondedNeighborhood(list &members, int traver // this->depthFirstSearch(members); } +// AS2023 - alternative call sig for logging that includes a log string +void Molecule::traverseBondedNeighborhood(list &members, int traversalLimit, string &logstr) +{ + //always call breadth first search, it is a bit faster + //if(traversalLimit>=0) + Molecule::breadthFirstSearch(members, this, traversalLimit, logstr); + //else + // this->depthFirstSearch(members); +} //Isn't ever called really, but is availabe. Note that it cannot use traversal limits //because it is depth first @@ -641,7 +801,6 @@ void Molecule::depthFirstSearch(list &members) } - void Molecule::printMoleculeList(list &members) { cout<<"List of molecules contains: "<defaultCompState = new int [numOfComponents]; this->isIntegerCompState = new bool [numOfComponents]; + int nostate = Molecule::NOSTATE; for(int c=0; ccompName[c]=compName.at(c); this->isIntegerCompState[c]=isIntegerComponent.at(c); @@ -543,6 +544,42 @@ void MoleculeType::updateRxnMembership(Molecule * m) rxn->tryToAdd(m, reactionPositions.at(r)); this->system->update_A_tot(rxn,oldA,rxn->update_a()); } + +} + +void MoleculeType::updateConnectedRxnMembership(Molecule * m, ReactionClass * firedReaction) +{ + // Replace the iteration over all reactions for the MoleculeType in + // MoleculeType::updateRxnMembership by only the + // connectedReactions for the fired Reaction. This is a much smaller loop + // and skips moleculetypes that are not the TemplateMolecule of the reactant + // in the connected reaction right away. + // Arvind Rasi Subramaniam + // + for (int r=0; rgetNumConnectedRxns(); r++) { + rxn = firedReaction->getconnectedRxn(r); + for (int pos=0; posgetNumOfReactants(); pos++) { + if (rxn->getMoleculeTypeOfReactantTemplate(pos) != this) continue; + double oldA = rxn->get_a(); + double oldAwithTotal = rxn->update_a(); + rxn->tryToAdd(m, pos); + this->system->update_A_tot(rxn,oldA,rxn->update_a()); + // Used for debugging to see which reaction rates changed + // upon updating molecule membership + // Arvind Rasi Subramaniam Nov 21, 2018 + if (!this->system->getTrackConnected()) continue; + double newA = rxn->update_a(); + if (oldAwithTotal != newA) { + this->system->getConnectedRxnFileStream() << + this->system->getGlobalEventCounter() << "\t" << + firedReaction->getName() << "\t" << + m->getMoleculeTypeName() << "\t" << + m->getUniqueID() << "\t" << + rxn->getName() << "\t" << + oldAwithTotal << "\t" << newA << endl; + } + } + } } @@ -755,16 +792,16 @@ void MoleculeType::printDetails() const // friend functions -template -NFstream& operator<<(NFstream& nfstream, const T& value) -{ - if (nfstream.useFile_) - nfstream.file_ << value; - else - nfstream.str_ << value; - - return nfstream; -} +// template +// NFstream& operator<<(NFstream& nfstream, const T& value) +// { +// if (nfstream.useFile_) +// nfstream.file_ << value; +// else +// nfstream.str_ << value; + +// return nfstream; +// } diff --git a/src/NFcore/observable.cpp b/src/NFcore/observable.cpp index 4e042afe..3765dbbd 100644 --- a/src/NFcore/observable.cpp +++ b/src/NFcore/observable.cpp @@ -135,8 +135,11 @@ void Observable::getTemplateMoleculeList(int &n_templates, TemplateMolecule **&t n_templates = this->n_templates; tmList = this->templateMolecules; } - - +// AS-2021 +void Observable::addReferenceToGlobalFunction(GlobalFunction *f) { + f->addCounterPointer(&count); +} +// AS-2021 void Observable::addReferenceToMyself(mu::Parser *p) { p->DefineVar(obsName,&count); diff --git a/src/NFcore/observable.hh b/src/NFcore/observable.hh index e1fc9062..a139b484 100644 --- a/src/NFcore/observable.hh +++ b/src/NFcore/observable.hh @@ -17,7 +17,9 @@ namespace NFcore class Molecule; class TemplateMolecule; class Complex; - + // AS-2021 + class GlobalFunction; + // AS-2021 //! Tracks the counts of predefined observables in the simulation. /*! @@ -99,6 +101,10 @@ namespace NFcore void addReferenceToMyself(mu::Parser *p); void addReferenceToMyself(string referenceName, mu::Parser *p); void addDependentRxn(ReactionClass *r); + + // AS-2021 + void addReferenceToGlobalFunction(GlobalFunction *f); + // AS-2021 virtual int isObservable(Molecule *m) const = 0; virtual int isObservable(Complex *c) const = 0; diff --git a/src/NFcore/reactionClass.cpp b/src/NFcore/reactionClass.cpp old mode 100644 new mode 100755 index 1dc29d84..a237b9a0 --- a/src/NFcore/reactionClass.cpp +++ b/src/NFcore/reactionClass.cpp @@ -13,7 +13,7 @@ ReactionClass::ReactionClass(string name, double baseRate, string baseRateParame { //cout<<"\n\ncreating reaction "<system=s; - this-> tagged = false; + this->tagged = false; totalRateFlag=false; isDimerStyle=false; @@ -30,7 +30,9 @@ ReactionClass::ReactionClass(string name, double baseRate, string baseRateParame //Set up the template molecules from the transformationSet this->n_reactants = transformationSet->getNreactants(); this->n_mappingsets = transformationSet->getNmappingSets(); - this->reactantTemplates = new TemplateMolecule *[n_reactants]; +// cout<<"n_reactants "<< this->n_reactants << endl; +// cout<<"n_mappingsets "<< this->n_mappingsets << endl; + this->reactantTemplates = new TemplateMolecule *[n_reactants]; vector tmList; vector hasMapGenerator; for(unsigned int r=0; rconnectedReactions.push_back(rxn); +} + +bool ReactionClass::isReactionConnected(ReactionClass * rxn) { + // First check if any of the operations share MoleculeType and components with + // one of the reactant templates of rxn. + return this->transformationSet->checkConnection(rxn); +} ReactionClass::~ReactionClass() { @@ -272,6 +283,21 @@ ReactionClass::~ReactionClass() delete [] mappingSet; delete [] isPopulationType; delete [] identicalPopCountCorrection; + connectedReactions.clear(); +} + +/** Fill the reactant and product templates for inferring reaction connectivity matrix + * @author Arvind Rasi Subramaniam + */ +void ReactionClass::setAllReactantAndProductTemplates(map reactants, + map products) { + map ::iterator it; + // Fill the reactant template pattern + for (it = reactants.begin(); it != reactants.end(); ++it) + this->allReactantTemplates.push_back(it->second); + // Fill the product template pattern + for (it = products.begin(); it != products.end(); ++it) + this->allProductTemplates.push_back(it->second); } @@ -304,23 +330,47 @@ void ReactionClass::resetBaseRateFromSystemParamter() { } +/** For use in MoleculeTye::updateRxnMembership + * @author Arvind Rasi Subramaniam + */ +MoleculeType *ReactionClass::getMoleculeTypeOfReactantTemplate(int pos) const { + // return reactantTemplates.at(pos)->getMoleculeType(); + return reactantTemplates[pos]->getMoleculeType(); +} + void ReactionClass::printDetails() const { cout<< name <<" (id="<rxnId<<", baseRate="<getReactantCount(r)<<" mappings|\t"; - cout<reactantTemplates[r]->getPatternString()<<"\n"; + n_mappings += this->getReactantCount(r); + } +// if (n_mappings == 0) return; + + cout << name << " (id=" << this->rxnId << ", baseRate=" << baseRate + << ", a=" << a << ", fired=" << fireCounter << " times )" << endl; + for (unsigned int r = 0; r < n_reactants; r++) { + cout << " -|" << this->getReactantCount(r) << " mappings|\t"; + cout << this->reactantTemplates[r]->getPatternString() << "\n"; //cout<<"head: "<reactantTemplates[r]->printDetails(cout); //reactantTemplates[r]->printDetails(); } - if(n_reactants==0) - cout<<" >No Reactants: so this rule either creates new species or does nothing."<No Reactants: so this rule either creates new species or does nothing." + << endl; + cout << "\n"; } - void ReactionClass::fire(double random_A_number) { + this->fire(random_A_number, false); +} + +// AS2023 - Alternative call signature to tell fire call when we are tracking +// each firing for the rxnlog argument +string ReactionClass::fire(double random_A_number, bool track) { //cout<FIRE "<checkMolecularity(mappingSet) ) { // wrong molecularity! this is a NULL event ++(System::NULL_EVENT_COUNTER); - return; + // AS2023 - we need to return a string now that this can return + // an event log if track is true + return string(""); } - // output something if the reaction was tagged - if(tagged) { - cout<<"#RT "<rxnId<<" "<system->getCurrentTime(); - for(unsigned int k=0; kgetNumOfMappings();p++) { - Molecule *mForTag = mappingSet[k]->get(p)->getMolecule(); - cout<<" "<getMoleculeTypeName()<getUniqueID(); - } - cout<<" ]"; - } - cout<getNumOfMappings();p++) { + // Molecule *mForTag = mappingSet[k]->get(p)->getMolecule(); + // cout<<" "<getMoleculeTypeName()<getUniqueID(); + // } + // cout<<" ]"; + // } + // cout<transformationSet->getListOfProducts(mappingSet,products,traversalLimit); - - // display product molecules for debugging.. - //for( molIter = products.begin(); molIter != products.end(); molIter++ ) { - // cout<<">>molecule: "<<(*molIter)->getMoleculeTypeName()<printDetails(); - // cout<<"<<"<onTheFlyObservables) { @@ -413,17 +454,21 @@ void ReactionClass::fire(double random_A_number) { } } - // Through the MappingSet, transform all the molecules as neccessary // This will also create new molecules, as required. As a side effect, // deleted molecules will be removed from observables. - this->transformationSet->transform(this->mappingSet); - + // AS2023 - if tracking is turned on, transform needs a string to build up + string logstr; + if (this->system->getReactionTrackingStatus()) { + logstr = this->transformationSet->transform(this->mappingSet, true); + + } else { + logstr = this->transformationSet->transform(this->mappingSet); + } // Add newly created molecules to the list of products this->transformationSet->getListOfAddedMolecules(mappingSet,products,traversalLimit); - // if complex bookkeeping is on, find all product complexes // (this is useful for updating Species Observables and TypeII functions, so keep the info handy). // NOTE: this is a brute force approach: check complex of each molecule. there may be a more @@ -471,7 +516,6 @@ void ReactionClass::fire(double random_A_number) { } } - // Now update reaction membership, functions, and update any DOR Groups // also, gather a list of typeII dependencies that will require updating typeII_products.clear(); @@ -491,10 +535,9 @@ void ReactionClass::fire(double random_A_number) { // NOTE: as a side-effect, DORreactions that depend on molecule-scoped local functions // (typeI relationship) will be updated as long as UTL is set appropriately. if ( mol->isAlive() ) - mol->updateRxnMembership(); + mol->updateRxnMembership(this, useConnectivity); } - // update complex-scoped local functions for typeII dependencies // NOTE: as a side-effect, dependent DOR reactions (via typeI molecule dependencies) will be updated if (system->getEvaluateComplexScopedLocalFunctions()) { @@ -531,27 +574,77 @@ void ReactionClass::fire(double random_A_number) { } } // done updating complex-scoped local functions - - // display final product molecules for debugging.. - //for( molIter = products.begin(); molIter != products.end(); molIter++ ) { - // cout<<">>molecule: "<<(*molIter)->getMoleculeTypeName()<printDetails(); - // cout<<"<<"<system->setLastRxnTime(this->system->getCurrentTime()); + // output to a JSON if the reaction was tagged + if (this->system->getReactionTrackingStatus()) { + if (tagged && track) { + string track_str = ""; + int level = 6; // indentation level + this->system->current_cpu_time = ((double) (clock() - this->system->start) / (double) CLOCKS_PER_SEC); + // we need the correct number of commas + if (this->system->getGlobalEventCounter() != 1) { + track_str += ",\n"; + } + // open firing and write info + track_str += std::string(level,' ') + "{\n" + + std::string(level+2,' ') + "\"props\": ["; + if (this->system->getRxnNumberTrack()) { + track_str += string("\"") + to_string(rxnId) + "\","; + } else { + track_str += string("\"") + name + "\","; + } + track_str += to_string(this->system->getGlobalEventCounter()) + + "," + to_string(this->system->getCurrentTime()) + "],\n"; + + // write transformation log + track_str += logstr; + + // close firing + track_str += std::string(level,' ') + "}"; + //Tidy up + products.clear(); + productComplexes.clear(); + return track_str; + } + } //Tidy up products.clear(); productComplexes.clear(); + // AS2023 - returning empty, if we are here logging was off + return ""; } +void ReactionClass::identifyConnectedReactions() { + ReactionClass * rxn; + vector allReactions; + allReactions = system->getAllReactions(); + for (unsigned int r=0; r < allReactions.size(); r++) { + rxn = allReactions.at(r); + if (this->isReactionConnected(rxn)) this->appendConnectedRxn(rxn); + } +} +bool ReactionClass::areMoleculeTypeAndComponentPresent(MoleculeType * mt, int cIndex) { + TemplateMolecule * t2; + for (unsigned int i=0; iisMoleculeTypeAndComponentPresent(mt, cIndex)) return true; + } + return false; +} +bool ReactionClass::isTemplateCompatible(TemplateMolecule * t) { + TemplateMolecule * t2; + for (unsigned int i=0; iisTemplateCompatible(t2)) return true; + } - - - - - + return false; +} diff --git a/src/NFcore/system.cpp b/src/NFcore/system.cpp index c44bcae6..a56c7c12 100644 --- a/src/NFcore/system.cpp +++ b/src/NFcore/system.cpp @@ -40,6 +40,8 @@ System::System(string name) ds=0; selector = 0; csvFormat = false; + anyRxnTagged = false; + max_cpu_time = -1; } @@ -67,6 +69,8 @@ System::System(string name, bool useComplex) ds=0; selector = 0; csvFormat = false; + anyRxnTagged = false; + max_cpu_time = -1; } System::System(string name, bool useComplex, int globalMoleculeLimit) @@ -92,6 +96,8 @@ System::System(string name, bool useComplex, int globalMoleculeLimit) ds=0; selector = 0; csvFormat = false; + anyRxnTagged = false; + max_cpu_time = -1; } @@ -276,9 +282,112 @@ void System::registerOutputFileLocation(string filename) } +/* + * Print reaction info if -rlog flag is given + * Note reactions can be tagged in BioNetGen or PySB or from the command line using the rtag flag + * @author: Rasi Subramaniam + */ +void System::registerReactionFileLocation(string filename) +{ + if (reactionOutputFileStream.is_open()) { reactionOutputFileStream.close(); } + reactionOutputFileStream.open(filename.c_str()); + + if(!reactionOutputFileStream.is_open()) { + cerr<<"Error in System! cannot open output stream to file "<=this->allReactions.size() ) { cerr<<"!!! Error when trying to tag reaction with reaction ID "<printDetails(); //<getName()==rName) + { + return (*rxnIter); + } + } + return 0; + cerr<<"!!! warning !!! cannot find reaction type name '"<< rName << "' in System: '"<name<<"'"<getMoleculeByUid(uid, true); +} +// AS2023 - alternative call sig to turn off warnings if we want to +Molecule * System::getMoleculeByUid(int uid, bool warn) { for( molTypeIter = allMoleculeTypes.begin(); molTypeIter != allMoleculeTypes.end(); molTypeIter++ ) { @@ -390,7 +520,9 @@ Molecule * System::getMoleculeByUid(int uid) return (*molTypeIter)->getMolecule(m); } } - cerr<<"!!! warning !!! cannot find active molecule with unique ID '"<< uid << "' in System: '"<name<<"'"<name<<"'"<prepareForSimulation(this); //cout<<"here 2..."<prepareForSimulation(this); //cout<<"here 3..."<getNumOfReactants()]; + + rxnIndexMap[r] = new int[allReactions.at(r)->getNumOfReactants()]; allReactions.at(r)->setRxnId(r); } + // Infer connected reactions if asked to do so from command line + // Arvind Rasi Subramaniam + if (connectivityFlag) { + // resize connected reactions map and intialize to false + connectedReactions = vector > (allReactions.size(), + vector (allReactions.size(), false)); + for(unsigned int r=0; ridentifyConnectedReactions(); + if ((r + 1) % 10 == 0) { + cout << "Connectivity inferred for " << r + 1 << " reactions." + << endl; + } + // prepare the connected reaction map for quick lookup + for (int r2 = 0; r2 < allReactions.at(r)->getNumConnectedRxns(); + r2++) { + int rxn2_id = + allReactions.at(r)->getconnectedRxn(r2)->getRxnId(); + connectedReactions[r][rxn2_id] = true; + } + // print connected reactions if given the switch + if (!this->getPrintConnected()) continue; + for (int r2 = 0; r2 < allReactions.at(r)->getNumConnectedRxns(); + r2++) { + this->getConnectedRxnListFileStream() << r << "\t" + << allReactions.at(r)->getName() << "\t" // << r << "\t" + << allReactions.at(r)->getconnectedRxn(r2)->getRxnId() + << "\t" + << allReactions.at(r)->getconnectedRxn(r2)->getName() + << endl; + } + } + } + + //cout<<"here 4..."<prepareForSimulation(); + } + //cout<<"here 7..."<getReactionTrackingStatus()) { + // start the JSON and write some info about the simulation + this->getReactionFileStream() << + "{" << endl << + " \"simulation\": {" << endl << + " \"info\": {" << endl << + " \"name\": \"" << this->getName() << "\"," << endl << + " \"global_molecule_limit\": " << to_string(this->getGlobalMoleculeLimit()) << "," << endl << + // " \"obs_count\": \"" << to_string(this->getMolObsCount()) << "\"," << endl << + " \"number_of_molecule_types\": " << to_string(this->getNumOfMoleculeTypes()) << "," << endl << + " \"number_of_molecules\": " << to_string(this->getNumOfMolecules()) << endl << + " }," << endl << + " \"molecule_types\": [" << endl; + + // prepare all molecule types + for(unsigned int mt=0; mtgetReactionFileStream() << + " {" << endl << + " \"name\": \"" + allMoleculeTypes.at(mt)->getName() + "\"," << endl << + " \"typeID\": " + to_string(allMoleculeTypes.at(mt)->getTypeID()) + ",\n" << + " \"components\": ["; + //deal with components + for (unsigned int mtci=0; mtcigetNumOfComponents(); mtci++) { + // we enter every component and a list of states + this->getReactionFileStream() << "\"" << + allMoleculeTypes.at(mt)->getComponentName(mtci) << "\""; + // deal with commas + if (mtci!=allMoleculeTypes.at(mt)->getNumOfComponents()-1) { + this->getReactionFileStream() << ","; + } + } + //close components and start component states + this->getReactionFileStream() << "],\n \"componentStates\": ["; + //deal with components + vector < vector < string > > comp_states = allMoleculeTypes.at(mt)->getPossibleCompStates(); + for (unsigned int mtci=0; mtcigetNumOfComponents(); mtci++) { + // we enter every component and a list of states + this->getReactionFileStream() << "["; + for (unsigned int mtcpsi=0; mtcpsigetReactionFileStream() << "\"" << comp_states[mtci][mtcpsi]; + // deal with commas + if (mtcpsi!=comp_states[mtci].size()-1) { + this->getReactionFileStream() << "\","; + } else { + this->getReactionFileStream() << "\""; + } + } + // deal with commas + if (mtci!=allMoleculeTypes.at(mt)->getNumOfComponents()-1) { + this->getReactionFileStream() << "],"; + } else { + this->getReactionFileStream() << "]"; + } + } + //close component states + this->getReactionFileStream() << "]\n }"; + // deal with commas + if (mt!=allMoleculeTypes.size()-1) { + this->getReactionFileStream() << ",\n"; + } else { + this->getReactionFileStream() << "\n"; + } + } + // close molecule types + this->getReactionFileStream() << + " ]," << endl; + + this->getReactionFileStream() << this->getSpeciesLog(); + + // close initial state and open firings for later + this->getReactionFileStream() << + " \"firings\": [" << endl; + } } @@ -663,6 +908,10 @@ double System::sim(double duration, long int sampleTimes, bool verbose) double end_time = current_time+duration; tryToDump(); + // AS2023 - depending on the tracking status we'll need a log string to build + string logstr; + bool logged = false; + while(current_timeprintAllObservableCounts(current_time); @@ -692,18 +941,26 @@ double System::sim(double duration, long int sampleTimes, bool verbose) curSampleTime+=dSampleTime; } if(verbose) { - cout << "Sim time: " << (curSampleTime-dSampleTime); - cout << "\tCPU time (total): " << ((double)(clock() - start)/(double)CLOCKS_PER_SEC) << "s"; - cout << "\t events (step): " << stepIteration< 0 & current_cpu_time > max_cpu_time) { + cout << "Max CPU time (" << max_cpu_time << ") reached, quitting." << endl; + break; + } } //cout<<"delta_t: " <getCount()<<"/"<getCount()<= r2*a_tot double randElement = getNextRxn(); @@ -715,37 +972,55 @@ double System::sim(double duration, long int sampleTimes, bool verbose) //nextReaction->printFullDetails(); //cout<printAllReactions(); - //this->printAllObservableCounts(this->current_time); - //cout<<"\n"; //Increment time iteration++; stepIteration++; globalEventCounter++; current_time+=delta_t; - - //5: Fire Reaction! (takes care of updates to lists and observables) - nextReaction->fire(randElement); - //this->printAllObservableCounts(this->current_time); - //cout<<"\n---"<printAllReactions(); +// cout << "Current reaction: " << "\n"; +// nextReaction->printDetails(); +// this->getMoleculeType(2)->getMolecule(0)->printDetails(); + // AS2023 - if we are tracking events, this needs to be dealt with here + if (this->getReactionTrackingStatus()) { + // AS2023 - getting the log for the event + logstr += nextReaction->fire(randElement, true); + // AS2023 - only write if we have a positive value for + // buffer size in events + if (this->getLogBufferSize()>0) { + // AS2023 - write if we have enough events stored in buffer + if ( (globalEventCounter % this->getLogBufferSize()) == 0) { + this->getReactionFileStream() << logstr; + // AS2023 - empty out the buffer + logstr = ""; + logged = true; + } + } + + } else { + nextReaction->fire(randElement); + } tryToDump(); - // outputAllPropensities(current_time, nextReaction->getRxnId()); - - //cout<getCount()<<"/"<getCount()<getName()=="Rule7") { - // cout<getCount()<<"/"<getCount()<getAllComplexes().printAllComplexes(); } if(curSampleTime-dSampleTime<(end_time-0.5*dSampleTime)) { outputAllObservableCounts(curSampleTime,globalEventCounter); } - + // AS2023 - if we missed a firing log, write what we have + if (!logged) { + this->getReactionFileStream() << logstr; + logstr = ""; + } + // Write list of molecule_types and reactions along with reaction firing counts + // TODO: Make this optional! + if (this->outputMoleculeTypesFile) { + outputAllMoleculeTypes(); + } + if (this->outputRxnFiringCountsFile) { + outputAllRxnFiringCounts(); + } finish = clock(); time = (double(finish)-double(start))/CLOCKS_PER_SEC; @@ -756,6 +1031,15 @@ double System::sim(double duration, long int sampleTimes, bool verbose) cout<<" Null events: "<< System::NULL_EVENT_COUNTER; cout<<" ("<<(time)/((double)iteration-(double)System::NULL_EVENT_COUNTER)<<" CPU seconds/non-null event )"<< endl; + // AS2023 - if we were tracking reactions, we should close the + // JSON file. We close the firing array, then the simulation + // level and finally the top level + if (this->getReactionTrackingStatus()) { + this->getReactionFileStream() << + "\n ]" << endl << + " }" << endl << + "}" << endl; + } cout.unsetf(ios::scientific); return current_time; } @@ -888,6 +1172,7 @@ void System::outputAllObservableNames() for(int k=0; kgetCount()); outputFileStream.write((char *) &count, sizeof(double)); } - if(outputGlobalFunctionValues) + if(outputGlobalFunctionValues) { for( functionIter = globalFunctions.begin(); functionIter != globalFunctions.end(); functionIter++ ) { + // AS-2021 + if ((*functionIter)->fileFunc==true) { + (*functionIter)->fileUpdate(); + } + // AS-2021 count=FuncFactory::Eval((*functionIter)->p); outputFileStream.write((char *) &count, sizeof(double)); } - + } if(outputEventCounter) { count=eventCounter; outputFileStream.write((char *) &count, sizeof(double)); @@ -1014,16 +1306,23 @@ void System::outputAllObservableCounts(double cSampleTime, int eventCounter) } else { if(!csvFormat) { - outputFileStream<<" "<getCount()); + outputFileStream<<"\t"<<((double)(*obsIter)->getCount()); } - if(outputGlobalFunctionValues) - for( functionIter = globalFunctions.begin(); functionIter != globalFunctions.end(); functionIter++ ) + if(outputGlobalFunctionValues) { + for( functionIter = globalFunctions.begin(); functionIter != globalFunctions.end(); functionIter++ ) { + // AS-2021 + if ((*functionIter)->fileFunc==true) { + (*functionIter)->fileUpdate(); + } + // AS-2021 outputFileStream<<" "<p); + } + } if(outputEventCounter) { - outputFileStream<<" "<getCount()); } - if(outputGlobalFunctionValues) - for( functionIter = globalFunctions.begin(); functionIter != globalFunctions.end(); functionIter++ ) + if(outputGlobalFunctionValues) { + for( functionIter = globalFunctions.begin(); functionIter != globalFunctions.end(); functionIter++ ) { + // AS-2021 + if ((*functionIter)->fileFunc==true) { + (*functionIter)->fileUpdate(); + } + // AS-2021 outputFileStream<<", "<p); + } + } if(outputEventCounter) { outputFileStream<<", "<getName(); @@ -1075,9 +1381,16 @@ void System::printAllObservableCounts(double cSampleTime,int eventCounter) cout<getCount(); - if(outputGlobalFunctionValues) - for( functionIter = globalFunctions.begin(); functionIter != globalFunctions.end(); functionIter++ ) - cout<<"\t"<p); + if(outputGlobalFunctionValues) { + for( functionIter = globalFunctions.begin(); functionIter != globalFunctions.end(); functionIter++ ) { + // AS-2021 + if ((*functionIter)->fileFunc==true) { + (*functionIter)->fileUpdate(); + } + // AS-2021 + cout<<"\t"<p)<getTypeID() << "\t" << (*molTypeIter)->getName() << endl; + } + moleculeTypeFileStream << this->getLastRxnTime() << "\tlast_rxn_firing_time" << endl; + moleculeTypeFileStream << this->current_time << "\tsimulated_time" << endl; + moleculeTypeFileStream << this->current_cpu_time << "\tcpu_time" << endl; + moleculeTypeFileStream.close(); +} + +void System::outputAllRxnFiringCounts() { + for(rxnIter = allReactions.begin(); rxnIter != allReactions.end(); rxnIter++ ) + { + rxnListFileStream << + (*rxnIter)->getRxnId() << "\t" << + (*rxnIter)->getFireCounter() << "\t" << (*rxnIter)->getName() << endl; + } + rxnListFileStream.close(); +} // NETGEN moved to ComplexList /* @@ -1403,7 +1737,7 @@ void System::evaluateAllLocalFunctions() { //Evaluate all local functions on this complex for(unsigned int l=0; levaluateOn(mol,LocalFunction::SPECIES); + localFunctions.at(l)->evaluateOn(mol,LocalFunction::SPECIES); //cout<<" value of function: "< +// NFstream& operator<<(NFstream& nfstream, const T& value) +// { +// if (nfstream.useFile_) +// nfstream.file_ << value; +// else +// nfstream.str_ << value; -// friend functions -template -NFstream& operator<<(NFstream& nfstream, const T& value) -{ - if (nfstream.useFile_) - nfstream.file_ << value; - else - nfstream.str_ << value; - - return nfstream; -} +// return nfstream; +// } diff --git a/src/NFcore/templateMolecule.cpp b/src/NFcore/templateMolecule.cpp index 04329864..c7cc5600 100644 --- a/src/NFcore/templateMolecule.cpp +++ b/src/NFcore/templateMolecule.cpp @@ -25,25 +25,25 @@ TemplateMolecule::TemplateMolecule(MoleculeType * moleculeType){ //Start everything off with no constraints, then we will reinitialize //every time we add some type of constraint - this->n_emptyComps=0; this->emptyComps=new int[0]; - this->n_occupiedComps=0; this->occupiedComps=new int[0]; - this->n_compStateConstraint=0; this->compStateConstraint_Comp=new int[0]; this->compStateConstraint_Constraint=new int[0]; - this->n_compStateExclusion=0; this->compStateExclusion_Comp=new int[0]; this->compStateExclusion_Exclusion=new int[0]; - - this->n_bonds=0; this->bondComp=new int[0]; this->bondCompName=new string[0]; this->bondPartner=new TemplateMolecule * [0]; this->bondPartnerCompName=new string[0]; this->bondPartnerCompIndex=new int[0]; this->hasVisitedBond=new bool[0]; + // + this->n_emptyComps=0; + this->n_occupiedComps=0; + this->n_compStateConstraint=0; + this->n_compStateExclusion=0; + this->n_bonds=0; this->n_connectedTo=0; this->connectedTo=new TemplateMolecule*[n_connectedTo]; this->hasTraversedDownConnectedTo=new bool[n_connectedTo]; @@ -81,6 +81,8 @@ TemplateMolecule::TemplateMolecule(MoleculeType * moleculeType){ //finally, we have to register this template molecule with the molecule //type so that we can easily destroy them at the end. this->moleculeType->addTemplateMolecule(this); + + this->mappedTm = NULL; } @@ -372,6 +374,11 @@ void TemplateMolecule::printDetails(ostream &o) { o<<"\n Map Generators: "; o<getPatternString(); + + o<<"\n Transformed Pattern: "; + o<getPatternString(); o<getCompIndexFromName(thisBsiteName); + + // TODO: these don't seem defined in merge, above block likely + // needs to be uncommented back //Copy over original the information to the new arrays for(int k=0; kgetCompIndexFromName(thisBsiteName); newBondComp[n_bonds] = compIndex; newBondCompName[n_bonds] = thisBsiteName; newBondPartner[n_bonds] = t2; @@ -514,7 +524,6 @@ void TemplateMolecule::addBond(string thisBsiteName, } newHasVisitedBond[n_bonds] = false; - //Delete the duplicated information delete [] bondComp; delete [] bondCompName; @@ -1177,7 +1186,7 @@ bool TemplateMolecule::compare(Molecule *m, ReactantContainer *rc, MappingSet *m //on this molecule Molecule *potentialMatch=m2->getBondedMolecule(bondPartnerCompIndex[b]); int thisBond = m2->getBondedMoleculeBindingSiteIndex(bondPartnerCompIndex[b]); - if(potentialMatch==Molecule::NOBOND) { + if(potentialMatch==nullptr) { //cout<<"potential match site has no bond"<getBondedMolecule(symBondPartnerCompIndex[c]); - if(potentialMatch==Molecule::NOBOND) { + if(potentialMatch==nullptr) { continue; } else if(potentialMatch!=matchMolecule) { continue; @@ -1521,6 +1530,138 @@ bool TemplateMolecule::compare(Molecule *m, ReactantContainer *rc, MappingSet *m } +/** To match two template molecules + * I will closely follow the 'compare' function above for + * comparing a TemplateMolecule and a Molecule. + * @author Arvind Rasi Subramaniam + * @param tm - Template molecule to match against + * @return true or false + */ +bool TemplateMolecule::isTemplateCompatible(TemplateMolecule * tm) { + // Make sure the TMs are of the same type + if(tm->getMoleculeType() != getMoleculeType()) { + return false; + } + + + // If there is no overlap between the components of the two TemplateMolecules, + // then they are compatible + vector allComps; + vector allComps_tm; + for (int i;in_emptyComps;i++) allComps_tm.push_back(tm->emptyComps[i]); + for (int i;in_occupiedComps;i++) allComps_tm.push_back(tm->occupiedComps[i]); + for (int i;in_bonds;i++) allComps_tm.push_back(tm->bondComp[i]); + // + + // Check each component of one TM against all components of other TM + bool compOverlap = false; + for (int i : allComps) { + if (find(allComps_tm.begin(), allComps_tm.end(), i) != allComps_tm.end()) { + compOverlap = true; + break; + } + } + // Still no overlap, so return + if (compOverlap == false) return true; + + + + vector ::iterator it; + // Check that sites that are occupied in one TemplateMolecule + // are not specified to be empty in the other TemplateMolecule + for (int i=0; i < n_occupiedComps; ++i) { + for(int j=0; jn_emptyComps; j++) { + if (tm->emptyComps[j] == occupiedComps[i]) return false; + } + } + // Check that sites that are empty in one TemplateMolecule + // are not specified to be occupied in the other TemplateMolecule + for (int i=0; i < n_emptyComps; ++i) { + for(int j=0; jn_occupiedComps; j++) { + if (tm->occupiedComps[j] == emptyComps[i]) return false; + } + } + // Check that sites that are empty in one TemplateMolecule + // are not bonded in the other TemplateMolecule + for (int i=0; i < n_emptyComps; ++i) { + for(int j=0; jn_bonds; j++) { + if (tm->bondComp[j] == emptyComps[i]) return false; + } + } + + // If a state is constrained, make sure that it is either not + // constrained or excluded in the target tm, or if it is, + // that the constraint or exclusion are not incompatible between + // the two templatemolecules + for (int i=0; i < n_compStateConstraint; ++i) { + for(int j=0; jn_compStateConstraint; j++) { + if (tm->compStateConstraint_Comp[j] == compStateConstraint_Comp[i]) { + if (compStateConstraint_Constraint[i] != tm->compStateConstraint_Constraint[j]) { + return false; + } + } + } + for(int j=0; jn_compStateExclusion; j++) { + if (tm->compStateExclusion_Comp[j] == compStateConstraint_Comp[i]) { + if (compStateConstraint_Constraint[i] != tm->compStateExclusion_Exclusion[j]) { + return false; + } + } + } + } + + // Repeat the above but now for all excluded states in the + // current TemplateMolecule + for (int i=0; i < n_compStateExclusion; ++i) { + for(int j=0; jn_compStateConstraint; j++) { + if (tm->compStateConstraint_Comp[j] == compStateExclusion_Comp[i]) { + if (compStateExclusion_Exclusion[i] != tm->compStateConstraint_Constraint[j]) { + return false; + } + } + } + for(int j=0; jn_compStateExclusion; j++) { + if (tm->compStateExclusion_Comp[j] == compStateExclusion_Comp[i]) { + if (compStateExclusion_Exclusion[i] != tm->compStateExclusion_Exclusion[j]) { + return false; + } + } + } + } + + // Make sure that components are bonded to the same partner at the same componentindex if + // the bonds are specified in both TemplateMolecules + for(int b=0; bn_emptyComps; j++) { + if (tm->emptyComps[j] == bondComp[b]) { + return false; + } + } + + //Check if this component is among the bondComps in the target TemplateMolecule + for(int j=0; jn_bonds; j++) { + if (tm->bondComp[j] == bondComp[b]) { + // find the bond number in the target molecule that matches bond Component + + // If the two bonding partners are not same molecule type, then return false. + // Assuming here that the second Template has a bond partner since it is + // the comp is in the list of bondComps. + if (bondPartner[b]->getMoleculeType() != tm->bondPartner[j]->getMoleculeType()) return false; + + // If the two bonding partners are not bonded on the same site, then return false + if (bondPartnerCompIndex[b] != tm->bondPartnerCompIndex[j]) return false; + } + } + } + + // We got this far, so all matches above did not raise any red flag + return true; +} string addStateConstraint(string original, string compName, string newConstraint) @@ -2005,9 +2146,23 @@ bool TemplateMolecule::checkSymmetryAroundBond(TemplateMolecule *tm1, TemplateMo return true; } - - - +bool TemplateMolecule::isMoleculeTypeAndComponentPresent(MoleculeType * mt, int cIndex) { + if (this->getMoleculeType() != mt) return false; + + // First make a joint vector of components specified in the TemplateMolecule + // TODO: ensure this is what is supposed to be done + for(int i=0; iemptyComps[i] == cIndex) return true; + } + for(int i=0; ioccupiedComps[i] == cIndex) return true; + } + for(int i=0; ibondComp[i] == cIndex) return true; + } + + return false; +} diff --git a/src/NFcore/templateMolecule.hh b/src/NFcore/templateMolecule.hh index 0c95ce66..3bd5ecdb 100644 --- a/src/NFcore/templateMolecule.hh +++ b/src/NFcore/templateMolecule.hh @@ -108,6 +108,13 @@ namespace NFcore Molecule *mappedFrom, string mappedFromComponent); bool isSymMapValid(); + /* Function to test whether two Template molecules match or compatible with each other. + * This is useful for testing if the reactants and products of a fired reaction will affect + * the reactants for another reaction + * @author: Arvind Rasi Subramaniam + */ + bool isTemplateCompatible(TemplateMolecule* tm); + ////////////////////////////////////////////////////////////////////////////////////////////// //returns false if they are not symmetric, or true if they are @@ -126,7 +133,12 @@ namespace NFcore void printPattern(); void printPattern(ostream &o); + // To get and set mapped reactant or product templatemolecule + // Arvind Rasi Subramaniam + void setMappedPartner(TemplateMolecule * tm) {mappedTm = tm;}; + TemplateMolecule * getMappedPartner() {return mappedTm;}; + bool isMoleculeTypeAndComponentPresent(MoleculeType * mt, int cIndex); protected: @@ -149,7 +161,7 @@ namespace NFcore // Which of the unique components must be empty (no bonds) int n_emptyComps; - int *emptyComps; + int *emptyComps; // Which of the unique components must be occupied (bonded to something, something // that is not specified) @@ -212,6 +224,10 @@ namespace NFcore static vector ::iterator tmVecIter; static list ::iterator tmIter; + // For tracking the reactant or product that this TemplateMolecule is + // transformed into + TemplateMolecule * mappedTm; + }; } diff --git a/src/NFfunction/NFfunction.hh b/src/NFfunction/NFfunction.hh index 3b3b9d54..9b91923d 100644 --- a/src/NFfunction/NFfunction.hh +++ b/src/NFfunction/NFfunction.hh @@ -179,6 +179,18 @@ namespace NFcore { return varRefTypes[varRefIndex]; } + // AS-2021 + void fileUpdate(); + double getCounterValue(); + void loadParamFile(string filePath); + void enableFileDependency(string FilePath); + void setCtrName(string name); + void addCounterPointer(double *count); + // unhooking system timer option for now + // void addSystemPointer(System *s); + bool fileFunc; + // AS-2021 + /*! This is the actual Parser object that keeps track of the function and has references to all of its arguments. It is publicly visible, so be careful with it! Use this variable to evaluate the function. @@ -197,6 +209,18 @@ namespace NFcore { unsigned int n_params; string *paramNames; + + // AS-2021 + string ctrType; + string ctrName; + // unhooking system timer option for now + // System *sysPtr; + int currInd; + int dataLen; + double *counter; + vector > data; + string filePath; + // AS-2021 }; @@ -373,6 +397,16 @@ namespace NFcore { void addTypeIMoleculeDependency(MoleculeType *mt); + // AS-2021 + void fileUpdate(); + double getCounterValue(); + void loadParamFile(string filePath); + void enableFileDependency(string FilePath); + void setCtrName(string name); + void addFunctionPointer(GlobalFunction *f); + bool fileFunc; + // AS-2021 + protected: @@ -414,6 +448,19 @@ namespace NFcore { double * refLfValues; mu::Parser *p; + + // AS-2021 + string ctrType; + string ctrName; + // unhooking system timer option for now + // System *sysPtr; + GlobalFunction *funcPtr; + int currInd; + int dataLen; + double *counter; + vector > data; + string filePath; + // AS-2021 }; diff --git a/src/NFfunction/compositeFunction.cpp b/src/NFfunction/compositeFunction.cpp index a8ae1b0d..69fef2af 100644 --- a/src/NFfunction/compositeFunction.cpp +++ b/src/NFfunction/compositeFunction.cpp @@ -48,6 +48,10 @@ CompositeFunction::CompositeFunction(System *s, } p=0; + + // AS-2021 + this->fileFunc = false; + // AS-2021 } CompositeFunction::~CompositeFunction() { @@ -215,6 +219,7 @@ void CompositeFunction::finalizeInitialization(System *s) this->refLfRefNames = new string[n_refLfs]; this->refLfScopes = new int[n_refLfs]; this->refLfValues = new double[n_refLfs]; + for(unsigned int i=0; irefLfInds[i]=lfIndexValues.at(i); this->refLfRefNames[i]=lfReferenceName.at(i); @@ -353,7 +358,13 @@ void CompositeFunction::printDetails(System *s) { cout<<" = "<originalExpression<fileFunc==true) { + gfs[f]->fileUpdate(); + } + // AS-2021 gfValues[f]=FuncFactory::Eval(gfs[f]->p); cout<<" global function: "<fileFunc==true) { + this->fileUpdate(); + } + // AS-2021 cout<<" Function last evaluated to: "<getMoleculeTypeName() << endl; + if(n_lfs>0) { //cout<<"evaluating composite function with local dependencies."<getNiceName()<<" with scope: "<refLfValues[i] = this->lfs[refLfInds[i]]->getValue(molList[refLfScopes[i]],scope[refLfScopes[i]]); + try{ + this->refLfValues[i] = this->lfs[refLfInds[i]]->getValue(molList[refLfScopes[i]],scope[refLfScopes[i]]); + } + catch (LocalFunctionException &lfe){ + //the parameter that we sent in is incorrect + lfe.setIndex(i); + throw lfe; + } //cout<<"answer: "<refLfValues[i]<fileFunc==true) { + this->fileUpdate(); + } + // AS-2021 + return FuncFactory::Eval(p); + //evaluate this function +} +// AS-2021 +void CompositeFunction::loadParamFile(string filePath) +{ + // setup our vectors + vector time; + vector values; + // open file for reading + ifstream file(filePath.c_str()); + // Report if file doesn't exist + if(!file.good()){ + cout<<"Error preparing function "<name<<" in class GlobalFunction!!"<> a >> b) { + // convert a to double + istringstream aos(a); + double d; + aos >> d; + // add it to time + time.push_back(d); + // convert b to double + istringstream bos(b); + bos >> d; + // add it to values + values.push_back(d); + } + // put the vectors into data vector + this->data.push_back(time); + this->data.push_back(values); + } catch (exception const & e) { + cout<<"Error preparing function "<name<<" in class GlobalFunction!!"<ctrType = "Function"; + // this->setCtrName(fPtr->getName()); + this->setCtrName("__TFUN__VAL__"); + this->funcPtr = fPtr; +} - //evaluate this function - return FuncFactory::Eval(p); +void CompositeFunction::setCtrName(string name) { + this->ctrName = name; +} + +void CompositeFunction::enableFileDependency(string filePath) { + // load file + // cout<<"file dependency of function: "<loadParamFile(filePath); + } catch (exception const & e) { + cout<<"Error preparing function "<filePath = filePath; + // this sets it up so that this function knows it's supposed + // to be pulling values from a file + this->fileFunc = true; + // initialize internal index + this->currInd = 0; + // pull data lenght so we can reuse it + this->dataLen = data[0].size(); +} + +double CompositeFunction::getCounterValue() { + // depending on the type of the observable counter + // get the actual value + double ctrVal; + if (ctrType == "Function") { + ctrVal = FuncFactory::Eval(this->funcPtr->p); + } + // unhooking system timer option for now + // else { + // // not sure but this is likely slower + // ctrVal = this->sysPtr->getCurrentTime(); + // } + return ctrVal; +} +void CompositeFunction::fileUpdate() { + // TODO: Error checking and reporting + + // get counter val + double ctrVal = this->getCounterValue(); + + // basic step function implementation + // if we got past the last point, keep returning + // the last point + if (currInd>dataLen-1) { + currInd = dataLen-1; + p->DefineConst(ctrName,data[1][currInd]); + return; + } else if (currInd==dataLen-1) { + p->DefineConst(ctrName,data[1][currInd]); + return; + } + // a simple way to do interval locating + if (data[0][currInd] < data[0][currInd+1]) { + // next point is higher than the current point, we + // are waiting for the counter value to be higher + // than our current point + + // return 0 if we don't have data yet + if(data[0][0]>=ctrVal) { + // we haven't gotten to the point where + // we can get a value out, return 0 + // cout<<"not there yet, returning 0"<DefineConst(ctrName,0); + return; + } + // go up by one if the counter value got past + // the next value in the array + if (ctrVal>=data[0][currInd+1]) { + currInd += 1; + } + // note that this makes no sense if they are equal + // TODO: Raise error if they are equal. Better yet, parse + // it ahead of time and make sure that doesn't happen + } else { + // next point is lower than the current point, we + // are waiting for the counter value to be lower + // than our current point + + // return 0 if we don't have data yet + if(data[0][0]<=ctrVal) { + // we haven't gotten to the point where + // we can get a value out, return 0 + // cout<<"not there yet, returning 0"<DefineConst(ctrName,0); + return; + } + // go up by one if the counter value got past + // the next value in the array + if (ctrVal<=data[0][currInd+1]) { + currInd += 1; + } + } + // // return value from the value array + p->DefineConst(ctrName,data[1][currInd]); + return; } +// AS-2021 \ No newline at end of file diff --git a/src/NFfunction/function.cpp b/src/NFfunction/function.cpp index fa3b4e33..f9fc4420 100644 --- a/src/NFfunction/function.cpp +++ b/src/NFfunction/function.cpp @@ -37,6 +37,10 @@ GlobalFunction::GlobalFunction(string name, this->paramNames[i]=paramNames.at(i); } p=0; + + // AS-2021 + this->fileFunc = false; + // AS-2021 } @@ -141,10 +145,173 @@ void GlobalFunction::printDetails() // } - if(p!=0) + if(p!=0) { + // AS-2021 + if (this->fileFunc==true) { + this->fileUpdate(); + } + // AS-2021 cout<<" Function currently evaluates to: "< time; + vector values; + // open file for reading + ifstream file(filePath.c_str()); + // Report if file doesn't exist + if(!file.good()){ + cout<<"Error preparing function "<name<<" in class GlobalFunction!!"<> a >> b) { + // convert a to double + istringstream aos(a); + double d; + aos >> d; + // add it to time + time.push_back(d); + // convert b to double + istringstream bos(b); + bos >> d; + // add it to values + values.push_back(d); + } + // put the vectors into data vector + this->data.push_back(time); + this->data.push_back(values); + } catch (exception const & e) { + cout<<"Error preparing function "<name<<" in class GlobalFunction!!"<ctrType = "Observable"; + this->counter = counter; +} + +void GlobalFunction::setCtrName(string name) { + this->ctrName = name; +} + +// unhooking system timer option for now +// void GlobalFunction::addSystemPointer(System *s) { +// this->ctrType = "System"; +// this->sysPtr = s; +// } + +void GlobalFunction::enableFileDependency(string filePath) { + // load file + // TODO: Err out if this fails + try { + this->loadParamFile(filePath); + } catch (exception const & e) { + cout<<"Error preparing function "<filePath = filePath; + // this sets it up so that this function knows it's supposed + // to be pulling values from a file + this->fileFunc = true; + // initialize internal index + this->currInd = 0; + // pull data lenght so we can reuse it + this->dataLen = data[0].size(); +} +double GlobalFunction::getCounterValue() { + // depending on the type of the observable counter + // get the actual value + double ctrVal; + if (ctrType == "Observable") { + ctrVal = (*counter); + } + // unhooking system timer option for now + // else { + // // not sure but this is likely slower + // ctrVal = this->sysPtr->getCurrentTime(); + // } + return ctrVal; } +void GlobalFunction::fileUpdate() { + // TODO: Error checking and reporting + // get counter val + double ctrVal = this->getCounterValue(); + // basic step function implementation + // if we got past the last point, keep returning + // the last point + if (currInd>dataLen-1) { + currInd = dataLen-1; + p->DefineConst(ctrName,data[1][currInd]); + return; + } else if (currInd==dataLen-1) { + p->DefineConst(ctrName,data[1][currInd]); + return; + } + // a simple way to do interval locating + if (data[0][currInd] < data[0][currInd+1]) { + // next point is higher than the current point, we + // are waiting for the counter value to be higher + // than our current point + + // return 0 if we don't have data yet + if(data[0][0]>=ctrVal) { + // we haven't gotten to the point where + // we can get a value out, return 0 + // cout<<"not there yet, returning 0"<DefineConst(ctrName,0); + return; + } + // go up by one if the counter value got past + // the next value in the array + if (ctrVal>=data[0][currInd+1]) { + currInd += 1; + } + // note that this makes no sense if they are equal + // TODO: Raise error if they are equal. Better yet, parse + // it ahead of time and make sure that doesn't happen + } else { + // next point is lower than the current point, we + // are waiting for the counter value to be lower + // than our current point + + // return 0 if we don't have data yet + if(data[0][0]<=ctrVal) { + // we haven't gotten to the point where + // we can get a value out, return 0 + // cout<<"not there yet, returning 0"<DefineConst(ctrName,0); + return; + } + // go up by one if the counter value got past + // the next value in the array + if (ctrVal<=data[0][currInd+1]) { + currInd += 1; + } + } + // // return value from the value array + p->DefineConst(ctrName,data[1][currInd]); + return; +} +// AS-2021 void GlobalFunction::printDetails(System *s) { @@ -159,9 +326,16 @@ void GlobalFunction::printDetails(System *s) cout<<" "<getParameter(paramNames[i])<fileFunc==true) { + cout<<" Function relies on file: "<filePath<fileUpdate(); + } + // AS-2021 cout<<" Function currently evaluates to: "<getLocalFunctionValue(typeI_localFunctionIndex.at(ti)); } } + LocalFunctionException lfe; + lfe.setType1_Mol(&typeI_mol); + throw lfe; } else if(scope==LocalFunction::MOLECULE) { //cout<<"Molecule scope."<setEvaluateComplexScopedLocalFunctions(evaluateComplexScopedLocalFunctions); + + // set inferring and using reaction connectivity flag + s->useConnectivityFlag(connectivityFlag); + //Read the key lists needed for the simulation and make sure they exist... TiXmlElement *pListOfParameters = pModel->FirstChildElement("ListOfParameters"); if(!pListOfParameters) { cout<<"\tNo 'ListOfParameters' tag found. Quitting."; delete s; return NULL; } @@ -94,6 +99,16 @@ System * NFinput::initializeFromXML( //(we do not enforce that functions must exist... yet) if(!pListOfFunctions) { cout<<"\tNo 'ListOfParameters' tag found. Quitting."; delete s; return NULL; } TiXmlElement *pListOfMoleculeTypes = pListOfParameters->NextSiblingElement("ListOfMoleculeTypes"); if(!pListOfMoleculeTypes) { cout<<"\tNo 'ListOfMoleculeTypes' tag found. Quitting."; delete s; return NULL; } + // we need to quit if we have compartments + TiXmlElement *pListOfCompartments = pListOfParameters->NextSiblingElement("ListOfCompartments"); + if(pListOfCompartments) { + // check to see if we have compartments + TiXmlElement *pCompElement; + pCompElement = pListOfCompartments->FirstChildElement("compartment"); + if (pCompElement) { + cout<<"\tCompartments aren't supported in NFsim. Quitting."; delete s; return NULL; + } + } TiXmlElement *pListOfSpecies = pListOfMoleculeTypes->NextSiblingElement("ListOfSpecies"); if(!pListOfSpecies) { cout<<"\tNo 'ListOfSpecies' tag found. Quitting."; delete s; return NULL; } TiXmlElement *pListOfReactionRules = pListOfSpecies->NextSiblingElement("ListOfReactionRules"); @@ -127,13 +142,18 @@ System * NFinput::initializeFromXML( if(!verbose) cout<<"-"; else cout<<"\n\tReading list of Species..."<setSpeciesLog(logstr); if(!verbose) cout<<"-"; else cout<<"\n\tReading list of Observables..."< ¶meter, @@ -617,6 +639,11 @@ bool NFinput::initStartSpecies( vector::iterator snIter; + // AS2023 - vectors to keep track of what's going on + // during initialization of the system + vector operations; + vector mgids; + vector mids; //Loop through all the species TiXmlElement *pSpec; @@ -626,7 +653,8 @@ bool NFinput::initStartSpecies( string speciesName; if(!pSpec->Attribute("id")) { cerr<<"Species tag without a valid 'id' attribute. Quiting"<Attribute("id"); } @@ -637,7 +665,8 @@ bool NFinput::initStartSpecies( string specCount; if(!pSpec->Attribute("concentration")) { cerr<<"Species "<Attribute("concentration"); } @@ -663,7 +692,8 @@ bool NFinput::initStartSpecies( } catch (std::runtime_error &e1) { if(parameter.find(specCount)==parameter.end()) { cerr<<"Could not find parameter: "<second; } @@ -674,7 +704,8 @@ bool NFinput::initStartSpecies( //Make sure we didn't try to create a negative number of molecules if(specCountInteger<0) { cerr<<"I cannot, in good conscience, make a negative number ("<FirstChildElement("ListOfMolecules"); if(!pListOfMol) { cerr<<"Species "<Attribute("name") || ! pMol->Attribute("id")) { cerr<<"!!!Error. Invalid 'Molecule' tag found when creating species '"<Attribute("name"); molUid = pMol->Attribute("id"); @@ -751,7 +785,8 @@ bool NFinput::initStartSpecies( { cerr << "!!!Error. Found mixed population and agent molecule types when creating species '" << speciesName << "'. Quitting"< usedComponentNames; @@ -768,7 +803,8 @@ bool NFinput::initStartSpecies( string compId,compName,compBondCount; if(!pComp->Attribute("id") || !pComp->Attribute("name") || !pComp->Attribute("numberOfBonds")) { cerr<<"!!!Error. Invalid 'Component' tag found when creating '"<Attribute("id"); compName = pComp->Attribute("name"); @@ -814,13 +850,15 @@ bool NFinput::initStartSpecies( } if(!couldPlaceSymComp) { cout<<"Too many symmetric sites specified, when creating species: "<genDefaultMolecule(); + // AS2023 - storing what has been generated, we need both the ID of the + // molecule type as well as the global ID that's assigned to the instance + mids.push_back(mol->getMoleculeType()->getTypeID()); + mgids.push_back(mol->getUniqueID()); //for(int i=0; isetComponentState((*snIter), (int)stateValue.at(k)); //} + + // AS2023 - this is here to reduce the number of operations written + // note that the default molecule starts the component state at the + // 0th state so any time we manually set that right after generating + // a default one, that's a reduntant operation + if ((int)stateValue.at(k)!=0) { + operations.push_back("[\"StateChange\"," + + to_string(mol->getUniqueID()) + "," + + to_string(mol->getMoleculeType()->getCompIndexFromName((*snIter))) + "," + + to_string((int)stateValue.at(k)) + "]"); + } + } molecules.at(molecules.size()-1).push_back(mol); @@ -933,7 +988,8 @@ bool NFinput::initStartSpecies( { cerr << "!! Attempt to create illegal bond in population species: " << speciesName << ". Quitting." << endl; - return false; + // AS2023 - fails now return empty strings + return ""; } //First get the information on the bonds in the complex @@ -943,7 +999,8 @@ bool NFinput::initStartSpecies( string bondId, bSite1, bSite2; if(!pBond->Attribute("id") || !pBond->Attribute("site1") || !pBond->Attribute("site2")) { cerr<<"!! Invalid Bond tag for species: "<Attribute("id"); bSite1 = pBond->Attribute("site1"); @@ -962,11 +1019,20 @@ bool NFinput::initStartSpecies( for(int j=0;jgetUniqueID()) + "," + + to_string(molecules.at(bSiteMolIndex1).at(j)->getMoleculeType()->getCompIndexFromName(bSiteName1.c_str())) + "," + + to_string(molecules.at(bSiteMolIndex2).at(j)->getUniqueID()) + "," + + to_string(molecules.at(bSiteMolIndex2).at(j)->getMoleculeType()->getCompIndexFromName(bSiteName2.c_str())) + "]"); } } catch (exception& e) { cout<<"!!!!Invalid site value for bond: '"<0){ + // AS2023 - to add compression we make the full molecule type + // vector first + int molec_size = *max_element(mgids.begin(), mgids.end()); + molec_size+=1; + vector molec_vec; + for(unsigned int isi=0; isiprintAllMoleculeTypes(); + // s->printAllMoleculeTypes(); - //If we got here, then we are indeed successful - return true; + // AS2023 - If we got here, then we are indeed successful + // and we are returning the log + return logstr; } catch (...) { cerr<<"Caught some unknown error when creating Species."< reaction_name_id_map; + int reaction_count = 0; for ( pRxnRule = pListOfReactionRules->FirstChildElement("ReactionRule"); pRxnRule != 0; pRxnRule = pRxnRule->NextSiblingElement("ReactionRule")) { @@ -1083,6 +1200,8 @@ bool NFinput::initReactionRules( //First, read in the template molecules using these data structures // maps reactant pattern ids to TemplateMolecule pointers map reactants; + // maps product pattern ids to TemplateMolecule pointers + map products; // maps component ids to component objects map comps; // points to TemplateMolecules for Reactants @@ -1111,6 +1230,7 @@ bool NFinput::initReactionRules( TiXmlElement *pListOfMols = pReactant->FirstChildElement("ListOfMolecules"); if(pListOfMols) { + /* At this point, only the first reactant molecule is sent back as a template - rasi */ TemplateMolecule *tm = readPattern(pListOfMols, s, parameter, allowedStates, reactantName, reactants, comps, symMap, verbose, suggestedTraversalLimit); if(tm==NULL) return false; templates.push_back(tm); @@ -1121,10 +1241,7 @@ bool NFinput::initReactionRules( } } - //Outputting all the templates for debugging purposes - //map::iterator it; - // for ( it=reactants.begin() ; it != reactants.end(); it++ ) - // cout << (*it).first << " => " << (*it).second->getMoleculeType()->getName() << endl; + /////////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -1301,7 +1418,67 @@ bool NFinput::initReactionRules( return false; } } + + // Read in the Product Patterns for this rule + TiXmlElement *pListOfProductPatterns = pRxnRule->FirstChildElement("ListOfProductPatterns"); + if(!pListOfProductPatterns) { + cout<<"!!!!!!!!!!!!!!!!!!!!!!!! Warning:: ReactionRule "<FirstChildElement("ProductPattern"); pProduct != 0; pProduct = pProduct->NextSiblingElement("ProductPattern")) + { + const char *productName = pProduct->Attribute("id"); + if(!productName) { + cerr<<"Product tag in reaction "<FirstChildElement("ListOfMolecules"); + if(pListOfMols) { + /* At this point, only the first reactant molecule is sent back as a template - rasi */ + readTemplatePattern(pListOfMols, s, allowedStates, productName, products, comps, symMap, verbose); + } + else { + cerr<<"Product pattern "<FirstChildElement("Map"); + if(!pListOfMaps) { + cout<<"!!!!!!!!!!!!!!!!!!!!!!!! Warning:: ReactionRule "<FirstChildElement("MapItem"); pMap != 0; pMap = pMap->NextSiblingElement("MapItem")) + { + // TODO: these don't have to exist + if ( !pMap->Attribute("sourceID") | !pMap->Attribute("targetID") ) { + continue; + } + reactantId = pMap->Attribute("sourceID"); + productId = pMap->Attribute("targetID"); + if ((reactantId.size() == 0) || (productId.size() == 0)) { + cerr<<"Map in reaction "<second->setMappedPartner(productIt->second); + productIt->second->setMappedPartner(reactantIt->second); + } + } /////////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -1452,6 +1629,9 @@ bool NFinput::initReactionRules( //Even though we had to make sure both ends exist, and we really only need one transformation //we give both templates so we can check for symmetric unbinding + // Arvind Rasi Subramaniam: Modified the method below to create + // an extra empty transformation for the binding partner + // for inferring connectivity later on. if(!ts->addUnbindingTransform(c1->t, c1->symPermutationName, c2->t, c2->symPermutationName)) return false; if(verbose) { cout<<"\t\t\t***Identified unbinding of site: "+c1->t->getMoleculeTypeName()+"("+c1->symPermutationName + ")"; @@ -1739,6 +1919,7 @@ bool NFinput::initReactionRules( //ts->finalize(); ReactionClass *r = 0; bool totalRateFlag=false; + bool tagFlag=false; /////////////////////////////////////////////////////////////////////////////////////////////////////////////////// // Read in the rate law for this reaction @@ -1763,6 +1944,21 @@ bool NFinput::initReactionRules( } } + // if( !pRateLaw->Attribute("tag") ) { + // cerr<<"\n!!Error! This XML file was generated using an older version of BioNetGen that does not support the 'TotalRate' convention!"<Attribute("tag")); + // if(rf>0) tagFlag=true; + // if(verbose) cout<<"\t\t\t= "<Attribute("id") || !pRateLaw->Attribute("type")) { cerr<<"!!Error:: ReactionRule "<addReaction(r); + //Finally, add the completed rxn rule to the system only + //base rate is non-zero. + if (r->getBaseRate() > 0) { + s->addReaction(r); + reaction_name_id_map[rxnName] = reaction_count; + reaction_count++; + } else { + if (verbose) { + cout << "\n!! Warning !! Rate Law " << r->getName() << + " not simulated due to zero base rate!!\n" << endl; + } + } + // Add the reactant and product templates to the reaction class + r->setAllReactantAndProductTemplates(reactants, products); r->setTotalRateFlag(totalRateFlag); + // if (tagFlag) { + r->tag(); + s->turnOnTagRxnOutput(); + // } + // Use reaction connectivity flag + // Set to true if given on the command line + r->setConnectivityFlag(s->getConnectivityFlag()); comps.clear(); } @@ -2563,7 +2778,7 @@ TemplateMolecule *NFinput::readPattern( if(allowedStates.find(molName+"_"+compName+"_"+compStateValue)==allowedStates.end()) { cerr<<"You are trying to give a pattern of type '"<second; @@ -2591,7 +2806,7 @@ TemplateMolecule *NFinput::readPattern( } catch (std::runtime_error e) { cerr<<"I couldn't parse the numberOfBonds value when creating pattern: "<second; @@ -2634,7 +2849,7 @@ TemplateMolecule *NFinput::readPattern( if(!lookup(symC, compId, comps, symMap)) { cerr<<"Could not find the symmetric component when creating a component state, but there\n"; cerr<<"should have been one!! I don't know what to do, so I'll quit."<symPermutationName); stateValue.push_back(stateValueInt); @@ -2666,7 +2881,7 @@ TemplateMolecule *NFinput::readPattern( } catch (std::runtime_error e) { cerr<<"I couldn't parse the numberOfBonds value when creating pattern: "<Attribute("id") || !pBond->Attribute("site1") || !pBond->Attribute("site2")) { cerr<<"!! Invalid Bond tag for pattern: "<Attribute("id"); bSite1 = pBond->Attribute("site1"); @@ -2799,7 +3014,7 @@ TemplateMolecule *NFinput::readPattern( cerr<<"!!!!Invalid site value for bond: '"< &allowedStates, + string patternName, + map &templates, + map &comps, + map &symMap, + bool verbose) +{ + try { + + //A vector to hold each template molecule as we compose the entire template molecule + vector < TemplateMolecule * > tMolecules; + + //Maps that map binding site ids into a molecule location in the molecules vector + //and the name of the binding site + map bSiteSiteMapping; + map bSiteMolMapping; + + //vectors that keep track of the states and thier specified values as we create the templates + vector stateName; + vector stateValue; + + vector emptyBondSite; + vector occupiedBondSite; + + //An iterator + vector ::iterator strVecIter; + + // Now loop through the molecules in the list + TiXmlElement *pMol; + bool foundTrash = false; + for ( pMol = pListOfMol->FirstChildElement("Molecule"); pMol != 0; pMol = pMol->NextSiblingElement("Molecule")) + { + //First get the type of molecule and retrieve the moleculeType object from the system + string molName, molUid; + if(!pMol->Attribute("name") || ! pMol->Attribute("id")) { + cerr<<"!!!Error. Invalid 'Molecule' tag found when creating pattern '"<Attribute("name"); + molUid = pMol->Attribute("id"); + } + + //Generate an error for any null or trash molecule + if(molName=="Null" || molName=="NULL" || molName=="null") { + // cerr<<"\n\nError! You cannot include a 'null' molecule in any reactant or observable pattern!"<getMoleculeTypeByName(molName); + TemplateMolecule *tempmol = new TemplateMolecule(moltype); + if(verbose) cout<<"\t\t\t\tIncluding Molecule of type: "< (molUid,c)); + + //Loop through the components of the molecule in order to set state values + TiXmlElement *pListOfComp = pMol->FirstChildElement("ListOfComponents"); + if(pListOfComp) + { + TiXmlElement *pComp; + for ( pComp = pListOfComp->FirstChildElement("Component"); pComp != 0; pComp = pComp->NextSiblingElement("Component")) + { + //Get the basic components of this molecule + string compId, compName, compBondCount; + if(!pComp->Attribute("id") || !pComp->Attribute("name") || !pComp->Attribute("numberOfBonds")) { + cerr<<"!!!Error. Invalid 'Component' tag found when creating '"<Attribute("id"); + compName = pComp->Attribute("name"); + compBondCount = pComp->Attribute("numberOfBonds"); + } + + //Set up basic component info so we can go back to it if need be... + component c(tempmol,compName); + comps.insert(pair (compId,c)); + + ////////////////////////////////////////////////////// + ////////////////////////////////////////////////////// + // Handle equivalent components off reaction center differently + // it is off reaction center if 1) it is an eq component and 2) it is not in the symMap + if(symMap.find(compId)==symMap.end() && moltype->isEquivalentComponent(compName)) { + //cout<<"we should treat this as a symmetric constraint"<Attribute("state")) { + string compStateValue = pComp->Attribute("state"); + if(compStateValue!="*" && compStateValue!="?") { + if(allowedStates.find(molName+"_"+compName+"_"+compStateValue)==allowedStates.end()) { + cerr<<"You are trying to give a pattern of type '"<second; + } + } + } + + //Check it as a binding site... + int bondConstraint = TemplateMolecule::NO_CONSTRAINT; + if(pComp->Attribute("numberOfBonds")) { + string numOfBonds = pComp->Attribute("numberOfBonds"); + int numOfBondsInt = -1; + + //const int MUST_BE_OCCUPIED = -2; + //const int EITHER_WAY_WORKS = -3; + if(numOfBonds.compare("+")==0) { + bondConstraint = TemplateMolecule::OCCUPIED; + } else if(numOfBonds.compare("*")==0) { + bondConstraint = TemplateMolecule::NO_CONSTRAINT; + } else if(numOfBonds.compare("?")==0) { + bondConstraint = TemplateMolecule::NO_CONSTRAINT; + } else { + try { + numOfBondsInt = NFutil::convertToInt(numOfBonds); + } catch (std::runtime_error e) { + cerr<<"I couldn't parse the numberOfBonds value when creating pattern: "<addSymCompConstraint(compName,compId,bondConstraint,stateConstraint); + ////////////////////////////////////////////////////// + ////////////////////////////////////////////////////// + ////////////////////////////////////////////////////// + } else { + //cout<<"can be mapped as a symmetric component, or is not a symmetric component"<Attribute("state")) { + string compStateValue = pComp->Attribute("state"); + //Make sure the given state is allowed (we allow for wildcards...) + if(compStateValue!="*" && compStateValue!="?") { + if(allowedStates.find(molName+"_"+compName+"_"+compStateValue)==allowedStates.end()) { + cerr<<"You are trying to give a pattern of type '"<second; + //Make sure we catch symmetric components in the case of a state change... + component *symC; //cout<<"state value: "<< compId; + if(!lookup(symC, compId, comps, symMap)) { + cerr<<"Could not find the symmetric component when creating a component state, but there\n"; + cerr<<"should have been one!! I don't know what to do, so I'll quit."<symPermutationName); + stateValue.push_back(stateValueInt); + } + } + } + //Check it as a binding site... + if(pComp->Attribute("numberOfBonds")) { + string numOfBonds = pComp->Attribute("numberOfBonds"); + int numOfBondsInt = -1; + //Only try to parse this bond as a number if the number + //of bonds is not the '+' character. The '+' character implies + //that the site is occupied without explicitly specifying who it is + //bound to. Now also handles the wild card character (implying that + //it can be bound or unbound - it doesn't matter. + const int MUST_BE_OCCUPIED = -2; + const int EITHER_WAY_WORKS = -3; + if(numOfBonds.compare("+")==0) { + numOfBondsInt = MUST_BE_OCCUPIED; + } else if(numOfBonds.compare("*")==0) { + numOfBondsInt = EITHER_WAY_WORKS; + } else if(numOfBonds.compare("?")==0) { + numOfBondsInt = EITHER_WAY_WORKS; + } else { + try { + numOfBondsInt = NFutil::convertToInt(numOfBonds); + } catch (std::runtime_error e) { + cerr<<"I couldn't parse the numberOfBonds value when creating pattern: "<symPermutationName); + } else if(numOfBondsInt==EITHER_WAY_WORKS) { + //add nothing if either way works of course! There + //is no constraint (the two ways are either bonded or not bonded) + } else if(numOfBondsInt==0) { + emptyBondSite.push_back(symC->symPermutationName); + } else if (numOfBondsInt==1) { + bSiteSiteMapping[compId] = symC->symPermutationName; + bSiteMolMapping[compId] = tMolecules.size(); + } else { + cerr<<"I can only handle a site that has 0 or 1 bonds in pattern: "<addComponentConstraint(*strVecIter,(int)stateValue.at(k)); + } + for(strVecIter = emptyBondSite.begin(); strVecIter != emptyBondSite.end(); strVecIter++ ) + { + tempmol->addEmptyComponent(*strVecIter); + } + for(strVecIter = occupiedBondSite.begin(); strVecIter != occupiedBondSite.end(); strVecIter++ ) + { + tempmol->addBoundComponent(*strVecIter); + } + + + //tempmol->printDetails(); + + + //Update our data storage with the new template and empty out the things we don't need + templates.insert(pair (molUid,tempmol)); + tMolecules.push_back(tempmol); + stateName.clear(); + stateValue.clear(); + emptyBondSite.clear(); + occupiedBondSite.clear(); + } + + //Here is where we add the bonds to the template molecules in the pattern + TiXmlElement *pListOfBonds = pListOfMol->NextSiblingElement("ListOfBonds"); + if(pListOfBonds) + { + //First get the information on the bonds in the complex + TiXmlElement *pBond; + for ( pBond = pListOfBonds->FirstChildElement("Bond"); pBond != 0; pBond = pBond->NextSiblingElement("Bond")) + { + //First, grab the bond information that we need to set things up + string bondId, bSite1, bSite2; + if(!pBond->Attribute("id") || !pBond->Attribute("site1") || !pBond->Attribute("site2")) { + cerr<<"!! Invalid Bond tag for pattern: "<Attribute("id"); + bSite1 = pBond->Attribute("site1"); + bSite2 = pBond->Attribute("site2"); + } + + //if(verbose)cout<<"reading bond "<::const_iterator iter = bSiteSiteMapping.begin(); +// iter != bSiteSiteMapping.end(); ++iter ) +// cout << iter->first << '\t' << iter->second << '\n'; +// cout<<"bSiteMolMapping"<::const_iterator iter = bSiteSiteMapping.begin(); +// iter != bSiteSiteMapping.end(); ++iter ) +// cout << iter->first << '\t' << iter->second << '\n'; + + + + //First look up the info from the component maps + if( bSiteSiteMapping.find(bSite1)!=bSiteSiteMapping.end() && + bSiteMolMapping.find(bSite1)!=bSiteMolMapping.end() && + bSiteSiteMapping.find(bSite2)!=bSiteSiteMapping.end() && + bSiteMolMapping.find(bSite2)!=bSiteMolMapping.end() + ) { + + string bSiteName1 = bSiteSiteMapping.find(bSite1)->second; + int bSiteMolIndex1 = bSiteMolMapping.find(bSite1)->second; + string bSiteName2 = bSiteSiteMapping.find(bSite2)->second; + int bSiteMolIndex2 = bSiteMolMapping.find(bSite2)->second; + TemplateMolecule::bind(tMolecules.at(bSiteMolIndex1),bSiteName1.c_str(),bSite1, + tMolecules.at(bSiteMolIndex2),bSiteName2.c_str(),bSite2); + + //Erase the bonds to make sure we don't add them again + bSiteSiteMapping.erase(bSite1); + bSiteMolMapping.erase(bSite1); + bSiteSiteMapping.erase(bSite2); + bSiteMolMapping.erase(bSite2); + } else { + + cerr<<"here"<::iterator strMapit; + for(strMapit = bSiteSiteMapping.begin(); strMapit != bSiteSiteMapping.end(); strMapit++ ) { + string bSiteId = (*strMapit).first; + string bSiteName = (*strMapit).second; + int bSiteMolIndex = bSiteMolMapping.find(bSiteId)->second; + + //skip symmetric sites here, because we already considered them earlier... + if(tMolecules.at(bSiteMolIndex)->getMoleculeType()->isEquivalentComponent(bSiteName.c_str())) { + continue; + } + + tMolecules.at(bSiteMolIndex)->addBoundComponent(bSiteName); + } + } catch (exception& e) { + cerr<<"Something went wacky when I was parsing pattern '"<printDetails(cout); + //} + + + //Now we have to find disjointed sets - that is whenever we have a Template + //Molecule that is connected, but not explicitly through bonds, we have to + //connect them via the connectedTo specification + + + //cout<<"checking for disjoint sets..."< > sets; + vector uniqueSetId; + int setCount = TemplateMolecule::getNumDisjointSets(tMolecules,sets,uniqueSetId); + +// cout<<"Unique Set Ids for the templates: "<1) { + // Possibly, we might want to enforce complex bookkeeping for such reactions.... + //if(!s->isUsingComplex()) { + // cout.flush(); + // cerr<<"Disjoint pattern found, but complex bookkeeping is turned off!"<getN_connectedTo(); + int ctIndex2=tMolecules.at(tm2)->getN_connectedTo(); + tMolecules.at(tm1)->addConnectedTo(tMolecules.at(tm2),ctIndex2); + tMolecules.at(tm2)->addConnectedTo(tMolecules.at(tm1),ctIndex1); + } + + + //for(unsigned int i=0; iprintDetails(); + //} + + //cout<<"traversing... let's see if we got everyone:"< tmList; + //TemplateMolecule::traverse(tMolecules.at(1),tmList); + //for(unsigned int i=0; iprintDetails(); + //} + + + } + + //Grab the first template molecule from the list, and arbitrarily set this as the root + if(tMolecules.empty()){ + if(foundTrash) { + // if(verbose) cout<<"\t\t\t\tWarning: You have an add molecule rule, but only a Trash or Null pattern listed..."< ¶meter, @@ -196,6 +199,22 @@ namespace NFinput { bool verbose, int &suggestedTraversalLimit); + //! Reads a pattern XML block for identifying product patterns + /*! Introduced to get a list of all reactant and product patterns for each reaction + * + * @return 1 for no error and 0 with errors + * @author Arvind Rasi Subramaniam + */ + int readTemplatePattern( + TiXmlElement * pListOfMol, + System * s, + map &allowedStates, + string patternName, + map &templates, + map &comps, + map &symMap, + bool verbose); + //! Reads a pattern XML block and returns the set of new TemplateMolecule objects. /*! @author Michael Sneddon diff --git a/src/NFinput/commandLineParser.cpp b/src/NFinput/commandLineParser.cpp index 337d4908..5b69d44d 100644 --- a/src/NFinput/commandLineParser.cpp +++ b/src/NFinput/commandLineParser.cpp @@ -52,7 +52,8 @@ bool NFinput::parseArguments(int argc, const char *argv[], map &a a++; } - //cout<<"found: '"<Attribute("type")) { + string funcType = pFunction->Attribute("type"); + string ctrType; + string ctrName; + if (funcType == "TFUN") { + // + if (!pFunction->Attribute("file")) { + cerr<<"!!!Error: TFUN type function "<Attribute("ctrName")) { + cerr<<"!!!Error: Can't find counter name for TFUN function "<Attribute("ctrName"); + } + // check references find our counter type + if(refNamesSorted.size()>0) { + TiXmlElement *refCheck; + string refCheckName; + string refCheckType; + for ( refCheck = pListOfRefs->FirstChildElement("Reference"); refCheck != 0; refCheck = refCheck->NextSiblingElement("Reference")) { + refCheckName = refCheck->Attribute("name"); + if ( refCheckName == ctrName ) { + // we found our counter + refCheckType = refCheck->Attribute("type"); + if ( refCheckType != "Observable" && refCheckType != "Function") { + cerr<<"!!!Error: TFUN type function "<Attribute("file"); + // we ensured we have the right type of ref name/type earlier + if (ctrType=="Observable") { + // make function file dependent + GlobalFunction *f = system->getGlobalFunctionByName(funcName); + f->enableFileDependency(filePath); + system->getObservableByName(ctrName)->addReferenceToGlobalFunction(f); + // f->setCtrName(refNamesSorted[0]); + f->setCtrName("__TFUN__VAL__"); + } else if (ctrType=="Function") { + CompositeFunction *f = system->getCompositeFunctionByName(funcName); + f->enableFileDependency(filePath); + // GlobalFunction *cfPtr = system->getGlobalFunctionByName(refNamesSorted[0]); + f->addFunctionPointer(system->getGlobalFunctionByName(ctrName)); + } else { + cerr<<"!!!Error: TFUN type function "<getGlobalFunctionByName(funcName); + // f->enableFileDependency(filePath); + // f->addSystemPointer(system); + // f->setCtrName("__COUNTER__"); + // } + } + } + // AS-2021 + //And here we clear our arrays argNames.clear(); refNames.clear(); diff --git a/src/NFreactions/mappings/mappingSet.cpp b/src/NFreactions/mappings/mappingSet.cpp index 5c71ac55..a0875d32 100644 --- a/src/NFreactions/mappings/mappingSet.cpp +++ b/src/NFreactions/mappings/mappingSet.cpp @@ -28,8 +28,9 @@ MappingSet::MappingSet(unsigned int id, vector &transformatio // add mapping with component index -1 mappings[t] = new Mapping(transformations.at(t)->getType(), -1 ); } - else + else { mappings[t] = new Mapping(transformations.at(t)->getType(), transformations.at(t)->getComponentIndex() ); + } } } MappingSet::~MappingSet() diff --git a/src/NFreactions/reactantLists/reactantList.cpp b/src/NFreactions/reactantLists/reactantList.cpp index f9809dc8..7dd7dd5d 100644 --- a/src/NFreactions/reactantLists/reactantList.cpp +++ b/src/NFreactions/reactantLists/reactantList.cpp @@ -159,6 +159,18 @@ MappingSet * ReactantList::getMappingSet(unsigned int mappingSetId) const { } +MappingSet * ReactantList::getWriteableMappingSet(unsigned int mappingSetId) { + + //Get the mappingSet position + int pos = msPositionMap[mappingSetId]; + + //might want to do some error checking here to be safe, but oh well. + + //give the mapping set back + return mappingSets[pos]; + +} + void ReactantList::removeMappingSet(unsigned int mappingSetId) { //Make sure this mappingSet is not empty diff --git a/src/NFreactions/reactantLists/reactantList.hh b/src/NFreactions/reactantLists/reactantList.hh index 4572d06f..69894781 100644 --- a/src/NFreactions/reactantLists/reactantList.hh +++ b/src/NFreactions/reactantLists/reactantList.hh @@ -90,7 +90,16 @@ namespace NFcore void pickRandomFromPopulation(MappingSet *&ms); + /** + * Changed to let the returned mappingset be modified. + * Using this in the case when a molecule is attached to a reaction + * to prevent its reordering in the reactant list + * @author Arvind Rasi Subramaniam + * @param mappingSetId + * @return + */ virtual MappingSet * getMappingSet(unsigned int mappingSetId) const; + virtual MappingSet * getWriteableMappingSet(unsigned int mappingSetId); /*! diff --git a/src/NFreactions/reactions/DORreaction.cpp b/src/NFreactions/reactions/DORreaction.cpp index 4b75aaa4..bd3372d6 100644 --- a/src/NFreactions/reactions/DORreaction.cpp +++ b/src/NFreactions/reactions/DORreaction.cpp @@ -102,7 +102,7 @@ DORRxnClass::DORRxnClass( if((unsigned)transform->getType()==TransformationFactory::LOCAL_FUNCTION_REFERENCE) { LocalFunctionReference *lfr = static_cast(transform); if(lfr->getPointerName()==lfArgumentPointerNameList.at(i)) { -// cout<<"Found a match here!"<getFunctionScope()<pushNextAvailableMappingSet(); + comparisonResult = reactantTemplates[reactantPos]->compare(m,reactantTree,ms,false,&symmetricMappingSet); if(!comparisonResult) { @@ -300,6 +301,7 @@ bool DORRxnClass::tryToAdd(Molecule *m, unsigned int reactantPos) { else{ //JJT: new mapping and we are keeping it, so evaluate the function and confirm the push double localFunctionValue = this->evaluateLocalFunctions(*it); + if(DEBUG_MESSAGE)cout<<"local function value is: "<confirmPush((*it)->getId(),localFunctionValue); m->setRxnListMappingId(rxnIndex,(*it)->getId()); @@ -352,6 +354,9 @@ bool DORRxnClass::tryToAdd(Molecule *m, unsigned int reactantPos) { } else { if(DEBUG_MESSAGE)cout<<"wasn't in the tree, so trying to push and compare"<pushNextAvailableMappingSet(); + if(DEBUG_MESSAGE)cout<<"calling comparsion method"<printDetails(); + comparisonResult = reactantTemplates[reactantPos]->compare(m,reactantTree,ms,false,&symmetricMappingSet); if(!comparisonResult) { if(DEBUG_MESSAGE)cout<<"shouldn't be in the tree, so we pop"<size(); } +/* +JJT: this function is called if the default mappingset information is sending the wrong parameter to the local function when using a species +scope label. for now the solution is to try every molecule referenced by the mapping set. This may be inefficient but it will only be as long +as the length of a pattern defined by the user (usually 5-6 mt long) multiplied by the lenght of the observable referenced by the local function, +so O(nm) with n, m <~ 6. +*/ +double DORRxnClass::pickLocalFunctionParameter(MappingSet* ms, int index, vector * type1_Mol, int* reactantCounts) +{ + for(unsigned int r=0; rgetNumOfMappings(); r++) + { + for (auto it: *(type1_Mol)){ + if(it == ms->get(r)->getMolecule()->getMoleculeType()){ + try{ + this->argMappedMolecule[index] = ms->get(r)->getMolecule(); + return this->cf->evaluateOn(this->argMappedMolecule,this->argScope, reactantCounts, this->n_reactants); + + } + catch(LocalFunctionException &lfe){ + if(lfe.getIndex() == index){ + continue; + } + else{ + return this->pickLocalFunctionParameter(ms, lfe.getIndex(), lfe.getType1_Mol(), reactantCounts); + } + + } + } + } + } + // there is truly no possible mapping. User mistake prob, the error message could use somem improvement + cout<<"Internal error in LocalFunction::evaluateOn()! Trying to evaluate a function with unknown scope."<n_argMolecules; i++) { //cout<<"here."<fileFunc==true) { + gf->fileUpdate(); + } + // AS-2021 a=FuncFactory::Eval(gf->p); } else if(cf!=0) { int * reactantCounts = new int[this->n_reactants]; @@ -111,9 +116,14 @@ void FunctionalRxnClass::printDetails() const { string trate = "off"; if(this->totalRateFlag) trate = "on"; - if(gf!=0) + if(gf!=0) { + // AS-2021 + if (gf->fileFunc==true) { + gf->fileUpdate(); + } + // AS-2021 cout<<"ReactionClass: " << name <<" ( baseFunction="<getNiceName()<<"="<p)<<", a="<n_reactants]; for(unsigned int r=0; rconnectivityFlag = s->getConnectivityFlag(); } @@ -255,6 +267,25 @@ int BasicRxnClass::checkForEquality(Molecule *m, MappingSet* ms, int rxnIndex, R } +/** + * Updates a molecule's reaction membership and a reaction's reactant list + * + * If a molecule matches the TemplateMolecule of a reaction, then it gets added. + * Note that NFsim has a single arbitrary TemplateMolecule for each reactant + * complex for a reaction. So it assumes that you will always check the molecule + * that matches the designated TemplateMolecule of a reaction for updates. So if + * I modify the identification of reactant molecules in + * TransformationSet::getListOfProducts(), I need to make sure that every molecule + * that is part of a reaction whose propensity might change is included in the + * listOfProducts. + * + * This is the gateway function to TemplateMolecule::compare(). + * @author Arvind Rasi Subramaniam + * @param m - molecule that is being compared against the TemplateMolecule of a reaction + * @param reactantPos - the reactant number among the reaction's reacant complexes + * (note that every connected complex gets only one reactant number.) + * @return true if there are no errors. + */ bool BasicRxnClass::tryToAdd(Molecule *m, unsigned int reactantPos) { //First a bit of error checking, that you should skip unless we are debugging... @@ -286,28 +317,63 @@ bool BasicRxnClass::tryToAdd(Molecule *m, unsigned int reactantPos) } } - /* - //Here we get the standard update... - if(m->getRxnListMappingId(rxnIndex)>=0) //If we are in this reaction... - { - if(!reactantTemplates[reactantPos]->compare(m)) { - // cout<<"Removing molecule "<getUniqueID()<<" which was at mappingSet: "<getRxnListMappingId(rxnIndex)<removeMappingSet(m->getRxnListMappingId(rxnIndex)); - m->setRxnListMappingId(rxnIndex,Molecule::NOT_IN_RXN); - } - - } else { - //Try to map it! - ms = rl->pushNextAvailableMappingSet(); - if(!reactantTemplates[reactantPos]->compare(m,rl,ms)) { - //we must remove, if we did not match. This will also remove - //everything that was cloned off of the mapping set - rl->removeMappingSet(ms->getId()); - } else { - m->setRxnListMappingId(rxnIndex,ms->getId()); - } - } - */ + // Here we get the standard update... + // if(m->getRxnListMappingId(rxnIndex)>=0) + // { + // if(!reactantTemplates[reactantPos]->compare(m)) { + // // cout<<"Removing molecule "<getUniqueID()<<" which was at mappingSet: "<getRxnListMappingId(rxnIndex)<removeMappingSet(m->getRxnListMappingId(rxnIndex)); + // m->setRxnListMappingId(rxnIndex,Molecule::NOT_IN_RXN); + // } + // // case when the molecule is not in the reaction + // } else { + // // Get a clean mappingSet from the reactantList + // // typically from the end: see the code for pusNextAvailableMappingSet() + // ms = rl->pushNextAvailableMappingSet(); + // if(!reactantTemplates[reactantPos]->compare(m,rl,ms)) { + // //we must remove, if we did not match. This will also remove + // //everything that was cloned off of the mapping set + // rl->removeMappingSet(ms->getId()); + // } else { + // m->setRxnListMappingId(rxnIndex,ms->getId()); + // } + // } + + // Arvind Rasi Subramaniam: I am modifying this so that the mappingSet does + // not change its position if the molecule still matches with the template. + // This will prevent the simulation trajectory being dependent on whether + // reaction-molecules pairs with unchanged membership are checked or not. + // if (connectivityFlag) { + // if(m->getRxnListMappingId(rxnIndex)>=0) + // { + // // Insted of removing the mappingSet + // // and then getting a different clean one from the reactant list, + // // get the mapping set corresponding to the molecule + // ms = rl->getWriteableMappingSet(m->getRxnListMappingId(rxnIndex)); + // // and clear it before use. + // ms->clear(); + // // If the molecule matches again, + // // the mappingSet gets updated during the compare routine. + // // If the molecule does not match anymore, remove the mapping set + // // and remove the rxn form the molecule's reaction membership. + // if(!reactantTemplates[reactantPos]->compare(m,rl,ms)) { + // rl->removeMappingSet(ms->getId()); + // m->setRxnListMappingId(rxnIndex,Molecule::NOT_IN_RXN); + // } + // // case when the molecule is not in the reaction + // } else { + // // Get a clean mappingSet from the reactantList + // // typically from the end: see the code for pusNextAvailableMappingSet() + // ms = rl->pushNextAvailableMappingSet(); + // if(!reactantTemplates[reactantPos]->compare(m,rl,ms)) { + // //we must remove, if we did not match. This will also remove + // //everything that was cloned off of the mapping set + // rl->removeMappingSet(ms->getId()); + // } else { + // m->setRxnListMappingId(rxnIndex,ms->getId()); + // } + // } + // } //Here we get the standard update... while(m->getRxnListMappingId(rxnIndex)>=0) { @@ -316,7 +382,6 @@ bool BasicRxnClass::tryToAdd(Molecule *m, unsigned int reactantPos) //m->setRxnListMappingId(rxnIndex,Molecule::NOT_IN_RXN); } - //Try to map it! ms = rl->pushNextAvailableMappingSet(); symmetricMappingSet.clear(); diff --git a/src/NFreactions/reactions/reaction.hh b/src/NFreactions/reactions/reaction.hh index 9db99e8d..3358f865 100644 --- a/src/NFreactions/reactions/reaction.hh +++ b/src/NFreactions/reactions/reaction.hh @@ -34,6 +34,8 @@ namespace NFcore protected: virtual void pickMappingSets(double randNumber) const; + // AS-6/22 + bool connectivityFlag; ReactantList **reactantLists; @@ -115,6 +117,8 @@ namespace NFcore virtual void pickMappingSets(double randNumber) const; + virtual double pickLocalFunctionParameter(MappingSet *ms, int, vector *, int*); + ReactantList **reactantLists; ReactantTree *reactantTree; diff --git a/src/NFreactions/transformations/moleculeCreator.cpp b/src/NFreactions/transformations/moleculeCreator.cpp index 06065053..aa52aaa2 100644 --- a/src/NFreactions/transformations/moleculeCreator.cpp +++ b/src/NFreactions/transformations/moleculeCreator.cpp @@ -93,6 +93,49 @@ MoleculeCreator::create_molecule() return molecule_object; } +// AS2023 - alternate call sig to store a lof of the molecule +// creation event for later writing +Molecule * +MoleculeCreator::create_molecule(string &logstr) +{ + + if ( isPopulationType() ) + { + // increment population + molecule_object->incrementPopulation(); + } + else + { + //Create the molecule + molecule_object = molecule_type->genDefaultMolecule(); + if (!logstr.empty()) { + logstr += " [\"Add\"," + + to_string(molecule_object->getUniqueID()) + + "," + to_string(molecule_object->getMoleculeType()->getTypeID()) + + "],\n"; + } + + //Set the component state values correctly + for( comp_iter = component_states.begin(); comp_iter != component_states.end(); ++comp_iter ) + { + molecule_object->setComponentState( (*comp_iter).first, (*comp_iter).second ); + if (!logstr.empty()) { + logstr += " [\"StateChange\"," + + to_string(molecule_object->getUniqueID()) + + "," + to_string((*comp_iter).first) + + "," + to_string((*comp_iter).second) + + "],\n"; + } + } + + //Prep the molecule and enterinto the simulation + molecule_type->addMoleculeToRunningSystemButDontUpdate( molecule_object ); + } + + // return a pointer to the new molecule, yeah! + return molecule_object; +} + Molecule * MoleculeCreator::get_population_pointer() const { diff --git a/src/NFreactions/transformations/moleculeCreator.hh b/src/NFreactions/transformations/moleculeCreator.hh index 61103810..92d5115a 100644 --- a/src/NFreactions/transformations/moleculeCreator.hh +++ b/src/NFreactions/transformations/moleculeCreator.hh @@ -33,6 +33,8 @@ namespace NFcore // create a molecule and get pointer Molecule * create_molecule( ); + // AS2023 - alternate call sig to write and store a log of molecule creation + Molecule * create_molecule( string &logstr ); // get the template molecule TemplateMolecule * getTemplateMolecule () const { return template_molecule; }; diff --git a/src/NFreactions/transformations/speciesCreator.cpp b/src/NFreactions/transformations/speciesCreator.cpp index eefd543b..1bbb1af1 100644 --- a/src/NFreactions/transformations/speciesCreator.cpp +++ b/src/NFreactions/transformations/speciesCreator.cpp @@ -246,3 +246,60 @@ void SpeciesCreator::create() //cout<genDefaultMolecule(); + if (!logstr.empty()) { + logstr += " [\"AddSpecies\"," + + to_string(newMoleculeCreations[m]->getUniqueID()) + + "," + to_string(newMoleculeCreations[m]->getMoleculeType()->getTypeID()) + + "],\n"; + } + } + + //Set the component state values correctly + for(int c=0; csetComponentState(ndStateIndex[c],ndStateValue[c]); + if (!logstr.empty()) { + logstr += string(" [\"StateChange\",") + + to_string(newMoleculeCreations[ndStateMolecule[c]]->getUniqueID()) + + "," + to_string(ndStateIndex[c]) + "," + + to_string(ndStateValue[c]) + "],\n"; + } + } + + //Make the bonds + for(unsigned int b=0; bgetUniqueID()) + + "," + to_string(bSite1[b]) + + "," + to_string(m2->getUniqueID()) + + "," + to_string(bSite2[b]) + "],\n"; + } + } + + //Create all the molecules + //Molecule *m = mt->genDefaultMolecule(); + + //Set all the states etc... + + //Prep the molecules and enter them into the simulation + for(unsigned int m=0; maddMoleculeToRunningSystem(newMoleculeCreations[m]); + //cout<<"Created:"<printDetails(); + //cout<cIndex=-1; + this->tm = NULL; +} +EmptyTransform::EmptyTransform(int cIndex) : + Transformation(TransformationFactory::EMPTY) +{ + this->cIndex=cIndex; + this->tm = NULL; +} +EmptyTransform::EmptyTransform(int cIndex, TemplateMolecule * tm) : + Transformation(TransformationFactory::EMPTY) +{ + this->cIndex=cIndex; + this->tm = tm; +} /////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////// @@ -13,11 +30,30 @@ StateChangeTransform::StateChangeTransform(int cIndex, int newValue) : { this->cIndex = cIndex; this->newValue = newValue; + this->tm = NULL; +} +StateChangeTransform::StateChangeTransform(int cIndex, int newValue, TemplateMolecule * tm) : + Transformation(TransformationFactory::STATE_CHANGE) +{ + this->cIndex = cIndex; + this->newValue = newValue; + this->tm = tm; } void StateChangeTransform::apply(Mapping *m, MappingSet **ms) { m->getMolecule()->setComponentState(cIndex,newValue); } +// AS2023 - alternative call sig to store a log of the transform +void StateChangeTransform::apply(Mapping *m, MappingSet **ms, string &logstr) +{ + m->getMolecule()->setComponentState(cIndex,newValue); + if (!logstr.empty()) { + logstr += string(" [\"StateChange\",") + + to_string(m->getMolecule()->getUniqueID()) + + "," + to_string(cIndex) + "," + + to_string(newValue) + "],\n"; + } +} @@ -27,12 +63,32 @@ IncrementStateTransform::IncrementStateTransform(unsigned int cIndex) : Transformation(TransformationFactory::INCREMENT_STATE) { this->cIndex = cIndex; + this->tm = NULL; +} +IncrementStateTransform::IncrementStateTransform(unsigned int cIndex, TemplateMolecule * tm) : + Transformation(TransformationFactory::INCREMENT_STATE) +{ + this->cIndex = cIndex; + this->tm = tm; } void IncrementStateTransform::apply(Mapping *m, MappingSet **ms) { int oldValue = m->getMolecule()->getComponentState(cIndex); m->getMolecule()->setComponentState(cIndex,oldValue+1); } +// AS2023 - alternative call sig to store a log of the transform +void IncrementStateTransform::apply(Mapping *m, MappingSet **ms, string &logstr) +{ + int oldValue = m->getMolecule()->getComponentState(cIndex); + m->getMolecule()->setComponentState(cIndex,oldValue+1); + if (!logstr.empty()) { + logstr += " [\"IncrementState\"," + + to_string(m->getMolecule()->getUniqueID()) + + "," + to_string(cIndex) + + "," + to_string(oldValue+1) + + "],\n"; + } +} /////////////////////////////////////////////////////////////// @@ -41,12 +97,32 @@ DecrementStateTransform::DecrementStateTransform(unsigned int cIndex) : Transformation(TransformationFactory::DECREMENT_STATE) { this->cIndex = cIndex; + this->tm = NULL; +} +DecrementStateTransform::DecrementStateTransform(unsigned int cIndex, TemplateMolecule * tm) : + Transformation(TransformationFactory::DECREMENT_STATE) +{ + this->cIndex = cIndex; + this->tm = tm; } void DecrementStateTransform::apply(Mapping *m, MappingSet **ms) { int oldValue = m->getMolecule()->getComponentState(cIndex); m->getMolecule()->setComponentState(cIndex,oldValue-1); } +// AS2023 - alternative call sig to store a log of the transform +void DecrementStateTransform::apply(Mapping *m, MappingSet **ms, string &logstr) +{ + int oldValue = m->getMolecule()->getComponentState(cIndex); + m->getMolecule()->setComponentState(cIndex,oldValue-1); + if (!logstr.empty()) { + logstr += " [\"DecrementState\"," + + to_string(m->getMolecule()->getUniqueID()) + + "," + to_string(cIndex) + + "," + to_string(oldValue-1) + + "],\n"; + } +} @@ -58,6 +134,10 @@ NewMoleculeBindingTransform::NewMoleculeBindingTransform(int cIndex, int otherRe BindingTransform(cIndex, otherReactantIndex, otherMappingIndex) { } +NewMoleculeBindingTransform::NewMoleculeBindingTransform(int cIndex, int otherReactantIndex, int otherMappingIndex, TemplateMolecule * tm) : + BindingTransform(cIndex, otherReactantIndex, otherMappingIndex, tm) +{ } + bool NewMoleculeBindingTransform::checkForNullCondition(Mapping *m, MappingSet **ms) { // binding to a new molecule cannot break a null condition (because null conditions are thus far @@ -74,9 +154,24 @@ BindingTransform::BindingTransform(int cIndex, int otherReactantIndex, int other this->cIndex=cIndex; this->otherReactantIndex=otherReactantIndex; this->otherMappingIndex=otherMappingIndex; + this->tm = NULL; +} +BindingTransform::BindingTransform(int cIndex, int otherReactantIndex, int otherMappingIndex, TemplateMolecule * tm) : + Transformation(TransformationFactory::BINDING) +{ + this->cIndex=cIndex; + this->otherReactantIndex=otherReactantIndex; + this->otherMappingIndex=otherMappingIndex; + this->tm = tm; } void BindingTransform::apply(Mapping *m, MappingSet **ms) +{ + Mapping *m2 = ms[this->otherReactantIndex]->get(this->otherMappingIndex); + Molecule::bind(m->getMolecule(),m->getIndex(), m2->getMolecule(), m2->getIndex()); +} +// AS2023 - alternative call sig to store a log of the transform +void BindingTransform::apply(Mapping *m, MappingSet **ms, string &logstr) { //cout<<" cIndex: "<getMolecule(),m->getIndex(), m2->getMolecule(), m2->getIndex()); + if (!logstr.empty()) { + logstr += " [\"AddBond\"," + + to_string(m->getMolecule()->getUniqueID()) + + "," + to_string(m->getIndex()) + + "," + to_string(m2->getMolecule()->getUniqueID()) + + "," + to_string(m2->getIndex()) + "],\n"; + } } - bool BindingTransform::checkForNullCondition(Mapping *m, MappingSet **ms) { // the null condition in this case is if a molecule is trying to bind a site to the same site on itself! @@ -136,10 +237,35 @@ UnbindingTransform::UnbindingTransform(int cIndex) : Transformation(TransformationFactory::UNBINDING) { this->cIndex=cIndex; + this->tm = NULL; +} +UnbindingTransform::UnbindingTransform(int cIndex, TemplateMolecule * tm) : + Transformation(TransformationFactory::UNBINDING) +{ + this->cIndex=cIndex; + this->tm = tm; } void UnbindingTransform::apply(Mapping *m, MappingSet **ms) +{ + vector tpl = Molecule::unbind(m->getMolecule(),m->getIndex()); + int m2id = tpl[0]; + int c2id = tpl[1]; +} +// AS2023 - alternative call sig to store a log of the transform +void UnbindingTransform::apply(Mapping *m, MappingSet **ms, string &logstr) { //cout<<"unbinding.."<getMolecule(),m->getIndex()); + vector tpl = Molecule::unbind(m->getMolecule(),m->getIndex()); + int m2id = tpl[0]; + int c2id = tpl[1]; + if (!logstr.empty()) { + logstr += " [\"DeleteBond\"," + + to_string(m->getMolecule()->getUniqueID()) + + "," + to_string(m->getIndex()) + + "," + to_string(m2id) + + "," + to_string(c2id) + + "],\n"; + } + } /////////////////////////////////////////////////////////////// @@ -148,19 +274,30 @@ AddSpeciesTransform::AddSpeciesTransform(SpeciesCreator *sc) : Transformation(TransformationFactory::ADD) { this->sc=sc; + this->tm = NULL; } +AddSpeciesTransform::AddSpeciesTransform(SpeciesCreator *sc, TemplateMolecule * tm) : + Transformation(TransformationFactory::ADD) +{ + this->sc=sc; + this->tm = tm; +} AddSpeciesTransform::~AddSpeciesTransform() { delete sc; } - void AddSpeciesTransform::apply(Mapping *m, MappingSet **ms) { this->sc->create(); } +// AS2023 - alternative call sig to store a log of the transform +void AddSpeciesTransform::apply(Mapping *m, MappingSet **ms, string &logstr) +{ + this->sc->create(logstr); +} /////////////////////////////////////////////////////////////// @@ -170,6 +307,15 @@ AddMoleculeTransform::AddMoleculeTransform( MoleculeCreator * _mc ) { this->mc = _mc; new_molecule = NULL; + this->tm = NULL; +} + +AddMoleculeTransform::AddMoleculeTransform( MoleculeCreator * _mc, TemplateMolecule * tm) + : Transformation( TransformationFactory::ADD ) +{ + this->mc = _mc; + new_molecule = NULL; + this->tm = tm; } @@ -193,7 +339,6 @@ AddMoleculeTransform::get_population_pointer() const return mc->get_population_pointer(); }; - void AddMoleculeTransform::apply_and_map(MappingSet *ms) { // create molecule and get pointer @@ -206,6 +351,19 @@ void AddMoleculeTransform::apply_and_map(MappingSet *ms) ms->set( im, new_molecule ); } } +// AS2023 - alternative call sig to store a log of the transform +void AddMoleculeTransform::apply_and_map(MappingSet *ms, string &logstr) +{ + // create molecule and get pointer + new_molecule = this->mc->create_molecule(logstr); + + // point mappings to the new molecule + unsigned int n_mappings = ms->getNumOfMappings(); + for ( unsigned int im = 0; im < n_mappings; ++im ) + { + ms->set( im, new_molecule ); + } +} @@ -214,14 +372,26 @@ void AddMoleculeTransform::apply_and_map(MappingSet *ms) RemoveMoleculeTransform::RemoveMoleculeTransform(int removalType) : Transformation(TransformationFactory::REMOVE) { this->removalType=removalType; + this->tm = NULL; } +RemoveMoleculeTransform::RemoveMoleculeTransform(int removalType, TemplateMolecule * tm) : + Transformation(TransformationFactory::REMOVE) { + this->removalType=removalType; + this->tm = tm; +} void RemoveMoleculeTransform::apply(Mapping *m, MappingSet **ms) { cout << "!! Warning: calling apply from a RemoveMoleculeTransform!" << "!! This cannot be handled here! The TransformationSet object should handle this!" << endl; } +// AS2023 - alternative call sig to store a log of the transform +void RemoveMoleculeTransform::apply(Mapping *m, MappingSet **ms, string &logstr) +{ + cout << "!! Warning: calling apply from a RemoveMoleculeTransform!" + << "!! This cannot be handled here! The TransformationSet object should handle this!" << endl; +} /////////////////////////////////////////////////////////////// @@ -230,11 +400,27 @@ DecrementPopulationTransform::DecrementPopulationTransform() : Transformation(TransformationFactory::DECREMENT_POPULATION) { this->cIndex = -1; + this->tm = NULL; +} +DecrementPopulationTransform::DecrementPopulationTransform(TemplateMolecule * tm) : + Transformation(TransformationFactory::DECREMENT_POPULATION) +{ + this->cIndex = -1; + this->tm = tm; } void DecrementPopulationTransform::apply(Mapping *m, MappingSet **ms) { m->getMolecule()->decrementPopulation(); } +// AS2023 - alternative call sig to store a log of the transform +void DecrementPopulationTransform::apply(Mapping *m, MappingSet **ms, string &logstr) +{ + m->getMolecule()->decrementPopulation(); + if (!logstr.empty()) { + logstr += " [\"DecrementPopulation\"," + + to_string(m->getMolecule()->getUniqueID()) + "],\n"; + } +} ///////////////////////////////////////////////////////////////// @@ -269,12 +455,6 @@ NFcore::Transformation * TransformationFactory::genNewMoleculeBindingTransform1( { return new NewMoleculeBindingTransform(bSiteIndex, otherReactantIndex, otherMappingIndex); } - -//NFcore::Transformation * TransformationFactory::genBindingSeparateComplexTransform1(unsigned int bSiteIndex, unsigned int otherReactantIndex, unsigned int otherMappingIndex) -//{ -// return new BindingSeparateComplexTransform(bSiteIndex, otherReactantIndex, otherMappingIndex); -//} - NFcore::Transformation * TransformationFactory::genBindingTransform2(unsigned int bSiteIndex) { return new EmptyTransform(bSiteIndex); @@ -319,6 +499,57 @@ Transformation * TransformationFactory::genLocalFunctionReference(string Pointer +NFcore::Transformation * TransformationFactory::genStateChangeTransform(unsigned int stateIndex, int newStateValue, TemplateMolecule * tm) +{ + return new StateChangeTransform(stateIndex, newStateValue, tm); +} +NFcore::Transformation * TransformationFactory::genBindingTransform1(unsigned int bSiteIndex, unsigned int otherReactantIndex, unsigned int otherMappingIndex, TemplateMolecule * tm) +{ + return new BindingTransform(bSiteIndex, otherReactantIndex, otherMappingIndex, tm); +} +NFcore::Transformation * TransformationFactory::genNewMoleculeBindingTransform1(unsigned int bSiteIndex, unsigned int otherReactantIndex, unsigned int otherMappingIndex, TemplateMolecule * tm) +{ + return new NewMoleculeBindingTransform(bSiteIndex, otherReactantIndex, otherMappingIndex, tm); +} +NFcore::Transformation * TransformationFactory::genBindingTransform2(unsigned int bSiteIndex, TemplateMolecule * tm) +{ + return new EmptyTransform(bSiteIndex, tm); +} +NFcore::Transformation * TransformationFactory::genUnbindingTransform(unsigned int bSiteIndex, TemplateMolecule * tm) +{ + return new UnbindingTransform(bSiteIndex, tm); +} +NFcore::Transformation * TransformationFactory::genUnbindingTransform2(unsigned int bSiteIndex, TemplateMolecule * tm) +{ + return new EmptyTransform(bSiteIndex, tm); +} +NFcore::AddSpeciesTransform * TransformationFactory::genAddSpeciesTransform(SpeciesCreator *sc, TemplateMolecule * tm) +{ + return new AddSpeciesTransform(sc, tm); +} +NFcore::AddMoleculeTransform * TransformationFactory::genAddMoleculeTransform(MoleculeCreator *mc, TemplateMolecule * tm) +{ + return new AddMoleculeTransform(mc, tm); +} +NFcore::Transformation * TransformationFactory::genRemoveMoleculeTransform(int removalType, TemplateMolecule * tm) +{ + return new RemoveMoleculeTransform(removalType, tm); +} + + +NFcore::Transformation * TransformationFactory::genIncrementStateTransform(unsigned int cIndex, TemplateMolecule * tm) +{ + return new IncrementStateTransform(cIndex, tm); +} +NFcore::Transformation * TransformationFactory::genDecrementStateTransform(unsigned int cIndex, TemplateMolecule * tm) +{ + return new DecrementStateTransform(cIndex, tm); +} + +NFcore::Transformation * TransformationFactory::genDecrementPopulationTransform(TemplateMolecule * tm) +{ + return new DecrementPopulationTransform(tm); +} diff --git a/src/NFreactions/transformations/transformation.hh b/src/NFreactions/transformations/transformation.hh index 30b9eebf..3a04fb53 100644 --- a/src/NFreactions/transformations/transformation.hh +++ b/src/NFreactions/transformations/transformation.hh @@ -27,6 +27,7 @@ namespace NFcore @author Michael Sneddon */ static Transformation * genStateChangeTransform(unsigned int cIndex, int newValue); + static Transformation * genStateChangeTransform(unsigned int cIndex, int newValue, TemplateMolecule * tm); /*! Generates a binding transformation for one of the two binding sites in the binding Transform. You will @@ -34,6 +35,7 @@ namespace NFcore @author Michael Sneddon */ static Transformation * genBindingTransform1(unsigned int bSiteIndex, unsigned int otherReactantIndex, unsigned int otherMappingIndex); + static Transformation * genBindingTransform1(unsigned int bSiteIndex, unsigned int otherReactantIndex, unsigned int otherMappingIndex, TemplateMolecule * tm); /*! Generates a binding transformation for one of the two binding sites in the binding Transform, where at least @@ -42,6 +44,7 @@ namespace NFcore @author Michael Sneddon */ static Transformation * genNewMoleculeBindingTransform1(unsigned int bSiteIndex, unsigned int otherReactantIndex, unsigned int otherMappingIndex); + static Transformation * genNewMoleculeBindingTransform1(unsigned int bSiteIndex, unsigned int otherReactantIndex, unsigned int otherMappingIndex, TemplateMolecule * tm); // deprecated @@ -55,6 +58,7 @@ namespace NFcore @author Michael Sneddon */ static Transformation * genBindingTransform2(unsigned int bSiteIndex); + static Transformation * genBindingTransform2(unsigned int bSiteIndex, TemplateMolecule * tm); /*! Generates an unbinding transformation for a particular binding site. Only @@ -63,6 +67,10 @@ namespace NFcore @author Michael Sneddon */ static Transformation * genUnbindingTransform(unsigned int bSiteIndex); + static Transformation * genUnbindingTransform(unsigned int bSiteIndex, TemplateMolecule* tm); + // added by Arvind Rasi Subramaniam for inferring reaction connectivity + static Transformation * genUnbindingTransform2(unsigned int bSiteIndex, TemplateMolecule * tm); + /*! Generates an Add Molecule transformation. @@ -75,7 +83,9 @@ namespace NFcore call the method "create_and_map" which is specifc to AddMoleculeTransform. */ static AddSpeciesTransform * genAddSpeciesTransform( SpeciesCreator * sc ); + static AddSpeciesTransform * genAddSpeciesTransform( SpeciesCreator * sc, TemplateMolecule * tm); static AddMoleculeTransform * genAddMoleculeTransform( MoleculeCreator * mc ); + static AddMoleculeTransform * genAddMoleculeTransform( MoleculeCreator * mc, TemplateMolecule * tm); /*! @@ -85,6 +95,7 @@ namespace NFcore @author Michael Sneddon */ static Transformation * genRemoveMoleculeTransform(int removalType); + static Transformation * genRemoveMoleculeTransform(int removalType, TemplateMolecule * tm); /*! Generates an empty transformation. This is used in cases where there is a reactant that is not transformed in a reaction, but that still needs @@ -98,12 +109,14 @@ namespace NFcore @author Michael Sneddon */ static Transformation * genIncrementStateTransform(unsigned int cIndex); + static Transformation * genIncrementStateTransform(unsigned int cIndex, TemplateMolecule * tm); /*! Generates an DecrementState transformation. @author Michael Sneddon */ static Transformation * genDecrementStateTransform(unsigned int cIndex); + static Transformation * genDecrementStateTransform(unsigned int cIndex, TemplateMolecule * tm); /*! Generates an IncrementPopulation transformation. @@ -117,6 +130,7 @@ namespace NFcore @author Justin Hogg */ static Transformation * genDecrementPopulationTransform(); + static Transformation * genDecrementPopulationTransform(TemplateMolecule * tm); /*! Indicates that a delete transform deletes the entire connected species */ @@ -185,13 +199,16 @@ namespace NFcore public: Transformation(int type) {this->type=type;}; virtual ~Transformation() {}; - int getType() const { return type; } + int getType() const { return type; }; virtual void apply(Mapping *m, MappingSet **ms) = 0; + // AS2023 - alternative call sig to store a log of the transform + virtual void apply(Mapping *m, MappingSet **ms, string & logstr) = 0; virtual int getComponentIndex() const = 0; virtual int getRemovalType() { return -1; }; // returns false if it does not meet a null condition, true if the reaction // should be rejected do to a null condition virtual bool checkForNullCondition(Mapping *m, MappingSet **ms) { return false; }; + virtual TemplateMolecule * getTemplateMolecule() const = 0; protected: int type; }; @@ -202,9 +219,12 @@ namespace NFcore LocalFunctionReference(string PointerName, int scope, TemplateMolecule *tm); virtual ~LocalFunctionReference() {}; virtual void apply(Mapping *m, MappingSet **ms) {}; + // AS2023 - alternative call sig to store a log of the transform + virtual void apply(Mapping *m, MappingSet **ms, string & logstr) {}; virtual int getComponentIndex() const { return -1; }; TemplateMolecule *getTemplateObject() const {return tm;}; + virtual TemplateMolecule * getTemplateMolecule() const {return this->tm;}; string getPointerName() const { return PointerName; }; int getFunctionScope() const {return scope; }; @@ -220,44 +240,60 @@ namespace NFcore class EmptyTransform : public Transformation { public: - EmptyTransform() : Transformation(TransformationFactory::EMPTY){ this->cIndex=-1; }; - EmptyTransform(int cIndex) : Transformation(TransformationFactory::EMPTY){ this->cIndex=cIndex; }; + EmptyTransform(); + EmptyTransform(int cIndex); + EmptyTransform(int cIndex, TemplateMolecule * tm); virtual ~EmptyTransform() {}; virtual void apply(Mapping *m, MappingSet **ms) {}; + // AS2023 - alternative call sig to store a log of the transform + virtual void apply(Mapping *m, MappingSet **ms, string & logstr) {}; virtual int getComponentIndex() const { return cIndex; }; + virtual TemplateMolecule * getTemplateMolecule() const {return this->tm;}; protected: int cIndex; + TemplateMolecule * tm; }; class StateChangeTransform : public Transformation { public: StateChangeTransform(int cIndex, int newValue); + StateChangeTransform(int cIndex, int newValue, TemplateMolecule * tm); virtual ~StateChangeTransform() {}; virtual void apply(Mapping *m, MappingSet **ms); + // AS2023 - alternative call sig to store a log of the transform + virtual void apply(Mapping *m, MappingSet **ms, string & logstr); virtual int getComponentIndex() const {return cIndex;}; + virtual TemplateMolecule * getTemplateMolecule() const {return this->tm;}; protected: int cIndex; int newValue; + TemplateMolecule * tm; }; class BindingTransform : public Transformation { public: BindingTransform(int cIndex, int otherReactantIndex, int otherMappingIndex); + BindingTransform(int cIndex, int otherReactantIndex, int otherMappingIndex, TemplateMolecule * tm); virtual ~BindingTransform() {}; virtual void apply(Mapping *m, MappingSet **ms); + // AS2023 - alternative call sig to store a log of the transform + virtual void apply(Mapping *m, MappingSet **ms, string & logstr); virtual int getComponentIndex() const {return cIndex;}; virtual bool checkForNullCondition(Mapping *m, MappingSet **ms); + virtual TemplateMolecule * getTemplateMolecule() const {return this->tm;}; protected: int cIndex; int otherReactantIndex; int otherMappingIndex; + TemplateMolecule * tm; }; class NewMoleculeBindingTransform : public BindingTransform { public: NewMoleculeBindingTransform(int cIndex, int otherReactantIndex, int otherMappingIndex); + NewMoleculeBindingTransform(int cIndex, int otherReactantIndex, int otherMappingIndex, TemplateMolecule * tm); virtual ~NewMoleculeBindingTransform() {}; virtual bool checkForNullCondition(Mapping *m, MappingSet **ms); }; @@ -275,23 +311,33 @@ namespace NFcore class UnbindingTransform : public Transformation { public: UnbindingTransform(int cIndex); + UnbindingTransform(int cIndex, TemplateMolecule * tm); virtual ~UnbindingTransform() {}; virtual void apply(Mapping *m, MappingSet **ms); + // AS2023 - alternative call sig to store a log of the transform + virtual void apply(Mapping *m, MappingSet **ms, string & logstr); virtual int getComponentIndex() const {return cIndex;}; + virtual TemplateMolecule * getTemplateMolecule() const {return this->tm;}; protected: int cIndex; + TemplateMolecule * tm; }; class AddSpeciesTransform : public Transformation { public: AddSpeciesTransform( SpeciesCreator * sc ); + AddSpeciesTransform( SpeciesCreator * sc , TemplateMolecule * tm); virtual ~AddSpeciesTransform(); - virtual void apply( Mapping *m, MappingSet **ms ); + virtual void apply( Mapping *m, MappingSet **ms); + // AS2023 - alternative call sig to store a log of the transform + virtual void apply( Mapping *m, MappingSet **ms, string & logstr); virtual int getComponentIndex() const {cerr<<"You should not get a component index from an AddMoleculeTransform!!"<tm;}; protected: SpeciesCreator * sc; + TemplateMolecule * tm; }; @@ -299,32 +345,44 @@ namespace NFcore { public: AddMoleculeTransform( MoleculeCreator * _mc ); + AddMoleculeTransform( MoleculeCreator * _mc , TemplateMolecule * tm); virtual ~AddMoleculeTransform(); virtual void apply( Mapping * m, MappingSet ** ms ) { cerr<<"apply method should not be called from an AddMoleculeTranform!!"<tm;}; protected: MoleculeCreator * mc; Molecule * new_molecule; + TemplateMolecule * tm; }; class RemoveMoleculeTransform : public Transformation { public: RemoveMoleculeTransform(int removalType); + RemoveMoleculeTransform(int removalType, TemplateMolecule * tm); virtual ~RemoveMoleculeTransform() {}; virtual void apply(Mapping *m, MappingSet **ms); + // AS2023 - alternative call sig to store a log of the transform + virtual void apply(Mapping *m, MappingSet **ms, string & logstr); virtual int getComponentIndex() const {cout<<"You should not get a component index from a RemoveMoleculeTransform!!"<tm;}; protected: int removalType; + TemplateMolecule * tm; }; @@ -333,32 +391,48 @@ namespace NFcore { public: DecrementPopulationTransform(); + DecrementPopulationTransform(TemplateMolecule * tm); virtual ~DecrementPopulationTransform() {}; virtual void apply(Mapping *m, MappingSet **ms); + // AS2023 - alternative call sig to store a log of the transform + virtual void apply(Mapping *m, MappingSet **ms, string & logstr); virtual int getComponentIndex() const { return cIndex; }; + virtual TemplateMolecule * getTemplateMolecule() const {return this->tm;}; + protected: int cIndex; + TemplateMolecule * tm; }; class IncrementStateTransform : public Transformation { public: IncrementStateTransform(unsigned int stateIndex); + IncrementStateTransform(unsigned int stateIndex, TemplateMolecule * tm); virtual ~IncrementStateTransform() {}; virtual void apply(Mapping *m, MappingSet **ms); + // AS2023 - alternative call sig to store a log of the transform + virtual void apply(Mapping *m, MappingSet **ms, string & logstr); virtual int getComponentIndex() const {return cIndex;}; + virtual TemplateMolecule * getTemplateMolecule() const {return this->tm;}; protected: int cIndex; + TemplateMolecule * tm; }; class DecrementStateTransform : public Transformation { public: DecrementStateTransform(unsigned int stateIndex); + DecrementStateTransform(unsigned int stateIndex, TemplateMolecule * tm); virtual ~DecrementStateTransform() {}; virtual void apply(Mapping *m, MappingSet **ms); + // AS2023 - alternative call sig to store a log of the transform + virtual void apply(Mapping *m, MappingSet **ms, string & logstr); virtual int getComponentIndex() const {return cIndex;}; + virtual TemplateMolecule * getTemplateMolecule() const {return this->tm;}; protected: int cIndex; + TemplateMolecule * tm; }; diff --git a/src/NFreactions/transformations/transformationSet.cpp b/src/NFreactions/transformations/transformationSet.cpp index fe2b4bdd..f7da1e14 100644 --- a/src/NFreactions/transformations/transformationSet.cpp +++ b/src/NFreactions/transformations/transformationSet.cpp @@ -1,16 +1,13 @@ - #include "transformationSet.hh" +#include "transformation.hh" using namespace NFcore; - - list TransformationSet::deleteList; list TransformationSet::updateAfterDeleteList; list ::iterator TransformationSet::it; - -TransformationSet::TransformationSet(vector reactantTemplates) +TransformationSet::TransformationSet(vector reactantTemplates) // @suppress("Class members should be properly initialized") { this->hasSymUnbinding=false; this->hasSymBinding = false; @@ -113,8 +110,7 @@ TransformationSet::~TransformationSet() } -TemplateMolecule * -TransformationSet::getTemplateMolecule( unsigned int reactantIndex ) const +TemplateMolecule * TransformationSet::getTemplateMolecule( unsigned int reactantIndex ) const { if ( reactantIndex < n_reactants ) { @@ -144,7 +140,7 @@ bool TransformationSet::addStateChangeTransform(TemplateMolecule *t, string cNam // 2) Create a Transformation object to remember the information //cout<<"Adding state change transform to value: "<getMoleculeType()->getCompIndexFromName(cName); - Transformation *transformation = TransformationFactory::genStateChangeTransform(cIndex, finalStateValue); + Transformation *transformation = TransformationFactory::genStateChangeTransform(cIndex, finalStateValue, t); // 3) Add the transformation object to the TransformationSet transformations[reactantIndex].push_back(transformation); @@ -190,7 +186,7 @@ bool TransformationSet::addIncrementStateTransform(TemplateMolecule *t, string c // 2) Create a Transformation object to remember the information int cIndex = t->getMoleculeType()->getCompIndexFromName(cName); - Transformation *transformation = TransformationFactory::genIncrementStateTransform(cIndex); + Transformation *transformation = TransformationFactory::genIncrementStateTransform(cIndex, t); // 3) Add the transformation object to the TransformationSet transformations[reactantIndex].push_back(transformation); @@ -215,7 +211,7 @@ bool TransformationSet::addDecrementStateTransform(TemplateMolecule *t, string c // 2) Create a Transformation object to remember the information int cIndex = t->getMoleculeType()->getCompIndexFromName(cName); - Transformation *transformation = TransformationFactory::genDecrementStateTransform(cIndex); + Transformation *transformation = TransformationFactory::genDecrementStateTransform(cIndex, t); // 3) Add the transformation object to the TransformationSet transformations[reactantIndex].push_back(transformation); @@ -266,11 +262,11 @@ bool TransformationSet::addBindingTransform(TemplateMolecule *t1, string bSiteNa //equal to the size. Transformation *transformation1; if(reactantIndex1==reactantIndex2) - transformation1 = TransformationFactory::genBindingTransform1(cIndex1, reactantIndex2, transformations[reactantIndex2].size()+1); + transformation1 = TransformationFactory::genBindingTransform1(cIndex1, reactantIndex2, transformations[reactantIndex2].size()+1, t1); else - transformation1 = TransformationFactory::genBindingTransform1(cIndex1, reactantIndex2, transformations[reactantIndex2].size()); + transformation1 = TransformationFactory::genBindingTransform1(cIndex1, reactantIndex2, transformations[reactantIndex2].size(), t1); - Transformation *transformation2 = TransformationFactory::genBindingTransform2(cIndex2); + Transformation *transformation2 = TransformationFactory::genBindingTransform2(cIndex2, t2); transformations[reactantIndex1].push_back(transformation1); MapGenerator *mg1 = new MapGenerator(transformations[reactantIndex1].size()-1); @@ -315,11 +311,11 @@ bool TransformationSet::addNewMoleculeBindingTransform(TemplateMolecule *t1, str //equal to the size. Transformation *transformation1; if(reactantIndex1==reactantIndex2) - transformation1 = TransformationFactory::genNewMoleculeBindingTransform1(cIndex1, reactantIndex2, transformations[reactantIndex2].size()+1); + transformation1 = TransformationFactory::genNewMoleculeBindingTransform1(cIndex1, reactantIndex2, transformations[reactantIndex2].size()+1, t1); else - transformation1 = TransformationFactory::genNewMoleculeBindingTransform1(cIndex1, reactantIndex2, transformations[reactantIndex2].size()); + transformation1 = TransformationFactory::genNewMoleculeBindingTransform1(cIndex1, reactantIndex2, transformations[reactantIndex2].size(), t1); - Transformation *transformation2 = TransformationFactory::genBindingTransform2(cIndex2); + Transformation *transformation2 = TransformationFactory::genBindingTransform2(cIndex2, t2); transformations[reactantIndex1].push_back(transformation1); MapGenerator *mg1 = new MapGenerator(transformations[reactantIndex1].size()-1); @@ -425,15 +421,26 @@ bool TransformationSet::addUnbindingTransform(TemplateMolecule *t, string bSiteN // 2) Create a Transformation object to remember the information unsigned int cIndex = tToTransform->getMoleculeType()->getCompIndexFromName(bSiteName); - Transformation *transformation = TransformationFactory::genUnbindingTransform(cIndex); + Transformation *transformation = TransformationFactory::genUnbindingTransform(cIndex, t); // 3) Add the transformation object to the TransformationSet transformations[reactantIndex].push_back(transformation); - // 3) Create a MapGenerator object and add it to the templateMolecule + // 4) Create a MapGenerator object and add it to the templateMolecule MapGenerator *mg = new MapGenerator(transformations[reactantIndex].size()-1); tToTransform->addMapGenerator(mg); + // 4) Create an empty transformation for the binding partner so that + // connectivity can be inferred. Arvind Rasi Subramaniam + if (t2 != 0) { + unsigned int cIndex2 = t2->getMoleculeType()->getCompIndexFromName(bSiteName2); + Transformation *transformation2 = TransformationFactory::genUnbindingTransform2(cIndex2, t2); + transformations[reactantIndex].push_back(transformation2); + MapGenerator *mg = new MapGenerator(transformations[reactantIndex].size()-1); + t2->addMapGenerator(mg); + } + + return true; } @@ -455,7 +462,7 @@ bool TransformationSet::addDeleteMolecule(TemplateMolecule *t, int deletionType) cerr<<" A(b).B(a),( instead of, say, A(b!1).B(a!1) ) you will get this error."<getTemplateMolecule() ); // 3) Add the transformation object to the TransformationSet addMoleculeTransformations.push_back( transformation ); @@ -558,9 +565,13 @@ int TransformationSet::find(TemplateMolecule *t) } return findIndex; } - - -bool TransformationSet::transform(MappingSet **mappingSets) +// AS2023 - normal calls should have tracking off +string TransformationSet::transform(MappingSet **mappingSets) +{ + return this->transform(mappingSets, false); +} +// AS2023 - alternative call sig to store a log of the transform +string TransformationSet::transform(MappingSet **mappingSets, bool tracking) { if(!finalized) { cerr<<"TransformationSet cannot apply a transform if it is not finalized!"<0) { for(int i=0; iapply_and_map( mappingSets[n_reactants+i] ); + // AS2023 - since we are in the tracking call, track the application + addMoleculeTransformations.at(i)->apply_and_map( mappingSets[n_reactants+i], logstr ); } } @@ -584,7 +603,8 @@ bool TransformationSet::transform(MappingSet **mappingSets) size = addSpeciesTransformations.size(); if(size>0) { for(int i=0; iapply(NULL,NULL); + // AS2023 - since we are in the tracking call, track the application + addSpeciesTransformations.at(i)->apply(NULL, NULL, logstr); } } @@ -599,16 +619,24 @@ bool TransformationSet::transform(MappingSet **mappingSets) Molecule * mol = ms->get(t)->getMolecule(); if ( transformations[r].at(t)->getRemovalType()==(int)TransformationFactory::COMPLETE_SPECIES_REMOVAL ) { // complex deletion: flag connected molecules for deletion - mol->traverseBondedNeighborhood(deleteList,ReactionClass::NO_LIMIT); + // AS2023 - since we are in the tracking call, track the deletion events + mol->traverseBondedNeighborhood(deleteList,ReactionClass::NO_LIMIT, logstr); } else { // molecule deletion: flag this molecule for deletion + // track deletions if tracking is on + // this has to be done here + if (tracking) { + // AS2023 - since we are in the tracking call, track the deletion operations + logstr += " [\"Delete\"," + to_string(mol->getUniqueID()) + "],\n"; + } deleteList.push_back( mol ); } } else { // handle other transforms - transformations[r].at(t)->apply(ms->get(t), mappingSets); + // AS2023 - since we are in the tracking call, track the operation + transformations[r].at(t)->apply(ms->get(t), mappingSets, logstr); } } } @@ -623,7 +651,13 @@ bool TransformationSet::transform(MappingSet **mappingSets) } deleteList.clear(); - return true; + // finalize tracking + if (tracking) { + // need to delete the last comma + logstr.erase(logstr.end()-2, logstr.end()); + logstr += "\n ]\n"; + } + return logstr; } @@ -677,7 +711,6 @@ bool TransformationSet::checkMolecularity( MappingSet ** mappingSets ) bool TransformationSet::getListOfProducts(MappingSet **mappingSets, list &products, int traversalLimit) { //if(!finalized) { cerr<<"TransformationSet cannot apply a transform if it is not finalized!"<::iterator molIter; for(unsigned int r=0; r &products, int traversalLimit) +// bool TransformationSet::getListOfAddedMolecules(MappingSet **mappingSets, vector &products, int traversalLimit) { //if(!finalized) { cerr<<"TransformationSet cannot apply a transform if it is not finalized!"<getTemplateMolecule(); + if (!t1) continue; + mt1 = t1->getMoleculeType(); + // AS2023 - if this is not a removal, track connections, removal + // doesn't give any reaction connections, so skip that + if (transfn->getType()!=(int)TransformationFactory::REMOVE) { + c1 = transfn->getComponentIndex(); + // If the moleculetype or component is not present in the other reaction, + // it is not connected + if (!rxn->areMoleculeTypeAndComponentPresent(mt1, c1)) continue; + + // If the TemplateMolecule is 'incompatible' with any of the reactants + // or products, then the reaction is not connected + if (!rxn->isTemplateCompatible(t1)) continue; + // Both checks passed for one op so return true + return true; + } else { + return false; + } + } + } + // Do the same as above but now for the transformed product templates + for(unsigned int r=0; rgetTemplateMolecule(); + if (!t1) continue; + t1 = t1->getMappedPartner(); + if (!t1) continue; + mt1 = t1->getMoleculeType(); + c1 = transfn->getComponentIndex(); + // If the moleculetype or component is present in the other reaction, + // it is not connected + if (!rxn->areMoleculeTypeAndComponentPresent(mt1, c1)) continue; + + // If the TemplateMolecule is 'incompatible' with any of the reactants + // or products, then the reaction is not connected + if (!rxn->isTemplateCompatible(t1)) continue; + // Both checks passed for one op so return true + return true; + } + } + // Now consider the newly added molecules where components don't matter, but + // compatibility + for (unsigned int i=0; igetTemplateMolecule(); + if (!t1) continue; + // If the TemplateMolecule is 'incompatible' with any of the reactants + // or products, then the reaction is not connected + if (!rxn->isTemplateCompatible(t1)) continue; + // Both checks passed for one op so return true + return true; + } + // Both checks did not pass for any reactant or product template, so not connected + return false; +} diff --git a/src/NFreactions/transformations/transformationSet.hh b/src/NFreactions/transformations/transformationSet.hh index 4bad0022..0bbb6967 100644 --- a/src/NFreactions/transformations/transformationSet.hh +++ b/src/NFreactions/transformations/transformationSet.hh @@ -17,6 +17,7 @@ namespace NFcore class Molecule; class SpeciesCreator; class MoleculeCreator; + class ReactionClass; //! Maintains a set of Transformation objects for a ReactionClass @@ -167,7 +168,9 @@ namespace NFcore MappingSet per reactant in the correct position in the array, please!). @author Michael Sneddon */ - bool transform(MappingSet **mappingSets); + string transform(MappingSet **mappingSets); + // AS2023 - alternative call sig when tracking is enabled + string transform(MappingSet **mappingSets, bool tracking); /*! Generates a blank MappingSet (blank in the sense that it is not mapped to @@ -297,6 +300,9 @@ namespace NFcore double getSymmetryFactor() const { return symmetryFactor; }; void setSymmetryFactor(double val) { symmetryFactor = val; useSymmetryFactor = true; }; + // To get the connected reactions for each transformation + bool checkConnection(ReactionClass * rxn); + protected: /*! diff --git a/src/NFscheduler/NFstream.cpp b/src/NFscheduler/NFstream.cpp index cc9d075c..6a105d1a 100644 --- a/src/NFscheduler/NFstream.cpp +++ b/src/NFscheduler/NFstream.cpp @@ -144,17 +144,17 @@ NFstream& NFstream::operator<<(NFstream& (*func)(NFstream &)) return ((*func))(*this); } -// friend functions -template -NFstream& operator<<(NFstream& nfstream, const T& value) -{ - if (nfstream.useFile_) - nfstream.file_ << value; - else - nfstream.str_ << value; - - return nfstream; -} +// // friend functions +// template +// NFstream& operator<<(NFstream& nfstream, const T& value) +// { +// if (nfstream.useFile_) +// nfstream.file_ << value; +// else +// nfstream.str_ << value; + +// return nfstream; +// } NFstream& endl(NFstream& nfstream) { diff --git a/src/NFscheduler/NFstream.h b/src/NFscheduler/NFstream.h index e38f6776..3ab40c28 100644 --- a/src/NFscheduler/NFstream.h +++ b/src/NFscheduler/NFstream.h @@ -53,8 +53,18 @@ class NFstream void check_mpi(); }; +// friend function definition in header to avoid separate instantiation template -NFstream& operator<< (NFstream& nfstream, const T& value); +NFstream& operator<< (NFstream& nfstream, const T& value) +{ + if (nfstream.useFile_) + nfstream.file_ << value; + else + nfstream.str_ << value; + + return nfstream; +} + NFstream& endl (NFstream& nfstream); diff --git a/src/NFsim.cpp b/src/NFsim.cpp index 76d490dd..7a194314 100644 --- a/src/NFsim.cpp +++ b/src/NFsim.cpp @@ -36,7 +36,7 @@ * number of features that make it ideal for handling large and complex * biochemical systems, such as functionally defined rate laws and reactions * that depend on local context. NFsim is designed to operate with the BioNetGen - * Language (http://bionetgen.org/). The new version of BNG is able to + * Language (http://bionetgen.org/). The new of BNG is able to * generate an XML encoded form of the BNG Language, which NFsim can take as input. * * For more details on setting up, running, and getting output from an NFsim simulation @@ -75,6 +75,38 @@ * -notf = disables On the Fly Observables, see manual * * -cb = turn on complex bookkeeping, see manual + * + * -connect - infer network connectivity before starting simulation. (default: no). + * @author Arvind Rasi Subramaniam + * + * -rxnlog [filename] - write out firing time and participating molecules for all reactions to a JSON file + * by default the expected extension is `.nfevent.json` + * @author Arvind Rasi Subramaniam + * + * -logbuffer [int] - how many firings to wait before writing to the rxnlog file + * Allows you to balance between CPU/memory impact of writing to a reaction log. + * @author Ali Sinan Saglam + * + * -trackconnected - write out the reactions whose rates change after firing of each reaction. + * this works only if -rxnlog switch is included + * @author: Arvind Rasi Subramaniam + * + * -printconnected - print connectivity of each reaction to an output file. (default: no). + * this works only if -rxnlog switch is included + * @author Arvind Rasi Subramaniam + * + * -trackrxnnum - track reaction number instead of name. this helps to keep the rxn log file small. + * this works only if -rxnlog switch is included + * @author: Arvind Rasi Subramaniam + * + * -maxcputime - maximum run time for simulation in seconds (default: no limit). + * @author Arvind Rasi Subramaniam + * + * -printmoltypes - output molecule types (default: false). + * @author Ali Sinan Saglam + * + * -printrxncounts - output reaction firing counts (default: false). + * @author Ali Sinan Saglam * * -gml [integer] = sets maximal number of molecules, per any MoleculeType, see manual * @@ -189,7 +221,7 @@ int main(int argc, char *argv[]) // turned off for the general release code. //if (!schedulerInterpreter(&argc, &argv)) return 0; - string versionNumber = "1.11"; + string versionNumber = "1.14.3"; cout<<"starting NFsim v"+versionNumber+"..."< argMap, bool verbose) globalMoleculeLimit = NFinput::parseAsInt(argMap,"gml",globalMoleculeLimit); } + bool connectivityFlag = false; + if (argMap.find("connect")!=argMap.end()) { + connectivityFlag = true; + } + //Actually create the system bool cb = false; if(turnOnComplexBookkeeping || blockSameComplexBinding) cb=true; int suggestedTraveralLimit = ReactionClass::NO_LIMIT; System *s = NFinput::initializeFromXML(filename,cb,globalMoleculeLimit,verbose, - suggestedTraveralLimit,evaluateComplexScopedLocalFunctions); + suggestedTraveralLimit, + evaluateComplexScopedLocalFunctions, + connectivityFlag); if(s!=NULL) @@ -466,6 +505,28 @@ System *initSystemFromFlags(map argMap, bool verbose) string outputFileName = argMap.find("o")->second; s->registerOutputFileLocation(outputFileName); s->outputAllObservableNames(); + if (argMap.find("printmoltypes")!=argMap.end()) { + s->registerMoleculeTypeFileLocation( + outputFileName.replace( + outputFileName.end()-5, + outputFileName.end(), + ".molecule_type_list.tsv")); + s->setOutputMoleculeTypes(true); + } else { + s->setOutputMoleculeTypes(false); + }; + + if (argMap.find("printrxncounts")!=argMap.end()) { + s->registerRxnListFileLocation( + outputFileName.replace( + outputFileName.end()-23, + outputFileName.end(), + ".rxn_list.tsv")); + s->setOutputRxnFiringCounts(true); + } else { + s->setOutputRxnFiringCounts(false); + }; + } else { if(s->isOutputtingBinary()) { s->registerOutputFileLocation(s->getName()+"_nf.dat"); @@ -475,8 +536,51 @@ System *initSystemFromFlags(map argMap, bool verbose) s->registerOutputFileLocation(s->getName()+"_nf.gdat"); s->outputAllObservableNames(); if(verbose) cout<<"\tStandard output will be written to: "<< s->getName()+"_nf.gdat" <registerMoleculeTypeFileLocation(s->getName() + "_molecule_type_list.tsv"); + s->registerRxnListFileLocation(s->getName() + "_rxn_list.tsv"); + } + } + + if (argMap.find("rxnlog") != argMap.end()) { + string rxnLogFileName = argMap.find("rxnlog")->second; + // AS2023 - register file location + s->registerReactionFileLocation(rxnLogFileName); + if (argMap.find("logbuffer") != argMap.end()) { + // AS2023 - set buffer size if given, 10k is the default + s->setLogBufferSize(stoul(argMap.find("logbuffer")->second)); + } + // track the reactions whose rates change upon each each reaction + // firing. This is useful for debugging to make sure that all the + // right reactions are updated after each firing. + // Arvind Rasi Subramaniam Nov 21, 2018 + if (argMap.find("trackconnected") != argMap.end()) { + s->registerConnectedRxnFileLocation( + rxnLogFileName.replace( + rxnLogFileName.end()-4, + rxnLogFileName.end(), + "_connected.tsv")); + s->setTrackConnected(true); + } else { + s->setTrackConnected(false); + } + if (argMap.find("printconnected") != argMap.end()) { + s->registerListOfConnectedRxnFileLocation( + rxnLogFileName.replace( + rxnLogFileName.end()-4, + rxnLogFileName.end(), + "_connectedlist.tsv")); + s->setPrintConnected(true); + } else { + s->setPrintConnected(false); + } + if (argMap.find("trackrxnnum") != argMap.end()) { + s->setRxnNumberTrack(true); + } else { + s->setRxnNumberTrack(false); } } + // } + //turn off on the fly calculation of observables if(argMap.find("notf")!=argMap.end()) { @@ -512,10 +616,17 @@ bool runFromArgs(System *s, map argMap, bool verbose) double eqTime = 0; double sTime = 10; int oSteps = 10; + double maxCpuTime = -1; //Get the simulation time that the user wants eqTime = NFinput::parseAsDouble(argMap,"eq",eqTime); sTime = NFinput::parseAsDouble(argMap,"sim",sTime); + + // if (argMap.find("maxcputime") != argMap.end()) { + // maxCpuTime = NFinput::parseAsDouble(argMap,"maxcputime",maxCpuTime); + // } + // s->setMaxCpuTime(maxCpuTime); + oSteps = NFinput::parseAsInt(argMap,"oSteps",(int)oSteps); //Prepare the system for simulation!! @@ -653,8 +764,43 @@ void printHelp(string version) cout<<" This allows you to run the same simulation and get the"< + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/test/Issue37/issue37_test2.bngl b/test/Issue37/issue37_test2.bngl new file mode 100644 index 00000000..940868a7 --- /dev/null +++ b/test/Issue37/issue37_test2.bngl @@ -0,0 +1,18 @@ +begin model + +begin species + A() 100 + B() 0 +end species + +begin observables + Molecules A A() +end observables + +begin reaction rules + 0 <-> A() 10,1 +end reaction rules + +end model + +simulate({method=>"nf",t_end=>100,n_steps=>100}) \ No newline at end of file diff --git a/test/Issue37/issue37_test2.xml b/test/Issue37/issue37_test2.xml new file mode 100644 index 00000000..4cf5dc48 --- /dev/null +++ b/test/Issue37/issue37_test2.xml @@ -0,0 +1,86 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/test/testSuite/t1.bngl b/test/testSuite/t1.bngl index 39629eac..8210d4d4 100644 --- a/test/testSuite/t1.bngl +++ b/test/testSuite/t1.bngl @@ -28,11 +28,11 @@ begin parameters end parameters -begin molecule type +begin molecule types CheY(p~unphos~PHOS) Motor(state~CW~CCW) -end molecule type +end molecule types begin seed species diff --git a/test/testSuite/t1.xml b/test/testSuite/t1.xml new file mode 100644 index 00000000..f3fc8966 --- /dev/null +++ b/test/testSuite/t1.xml @@ -0,0 +1,272 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + (MotCW*w)*exp(((g1/2.0)*((1.0/2.0)-((CheYp/((cellVolume*_Na)*(10^-6)))/(Kd+(CheYp/((cellVolume*_Na)*(10^-6)))))))) + + + + + + + + + + + + (MotCCW*w)*exp(((-(g1/2.0))*((1.0/2.0)-((CheYp/((cellVolume*_Na)*(10^-6)))/(Kd+(CheYp/((cellVolume*_Na)*(10^-6)))))))) + + + + diff --git a/test/testSuite/t2.bngl b/test/testSuite/t2.bngl index e6fe0218..9c1de9cf 100644 --- a/test/testSuite/t2.bngl +++ b/test/testSuite/t2.bngl @@ -40,9 +40,9 @@ begin observables end observables -begin function +begin functions #No functions in this model -end function +end functions diff --git a/test/testSuite/t3.bngl b/test/testSuite/t3.bngl index cea1ae4b..40ac5283 100644 --- a/test/testSuite/t3.bngl +++ b/test/testSuite/t3.bngl @@ -17,7 +17,7 @@ begin parameters kr 7 #per second rate of methylation kb 20 #per second rate of demethylation - kcat 10 #per second, base rate of CheA autophosphorylation + kCat 10 #per second, base rate of CheA autophosphorylation kDephos 10 #per second, rate of CheA dephosphorylation #initial number of receptors @@ -49,15 +49,15 @@ begin observables end observables -begin function +begin functions #Here we have just one local function, which just multiplies #the auto phosphorylation rate times the number of receptor sites #that are methylated. Note that observables are calculated on #the local context of parameter $1 - pOn($1) = kCat*(R1($1)+2*R2($1)+3*R3($1)+4*R4($1)) + pOn(x) = kCat*(R1(x)+2*R2(x)+3*R3(x)+4*R4(x)) -end function +end functions @@ -74,7 +74,7 @@ begin reaction rules #here is the actual local function based reaction. The local function pointer #points to the entire species (even though there is only one molecule in the species) - $1::Receptor(a!1).CheA(rec!1,p~unphos) -> $1::Receptor(a!1).CheA(rec!1,p~PHOS) pOn($1) + %x::Receptor(a!1).CheA(rec!1,p~unphos) -> %x::Receptor(a!1).CheA(rec!1,p~PHOS) pOn(x) #the dephos reaction of CheA is independent of the pOn function CheA(p~PHOS) -> CheA(p~unphos) kDephos diff --git a/test/testSuite/t3.xml b/test/testSuite/t3.xml new file mode 100644 index 00000000..b2cd95a9 --- /dev/null +++ b/test/testSuite/t3.xml @@ -0,0 +1,587 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + kCat*(((R1(x)+(2*R2(x)))+(3*R3(x)))+(4*R4(x))) + + + + + + + + + + pOn(x) + + + + diff --git a/test/testSuite/t6.bngl b/test/testSuite/t6.bngl index d90df4b4..d38dd33e 100644 --- a/test/testSuite/t6.bngl +++ b/test/testSuite/t6.bngl @@ -27,11 +27,11 @@ begin parameters end parameters -begin molecule type +begin molecule types CheY(p~unphos~PHOS) Motor(state~CW~CCW) -end molecule type +end molecule types begin seed species diff --git a/test/testSuite/t_dor2.bngl b/test/testSuite/t_dor2.bngl new file mode 100644 index 00000000..002e2205 --- /dev/null +++ b/test/testSuite/t_dor2.bngl @@ -0,0 +1,75 @@ +# t_dor2.bngl +# This model tests DOR2RxnClass functionality by creating reactions that +# depend on local functions with two different argument pointers. +# It models a bimolecular reaction between two different receptor types +# where the reaction rate depends on the methylation states of both reactants. + +begin parameters + + k_base 1.0 # base rate constant + + # initial receptor counts + RecA_count 100 + RecB_count 100 + +end parameters + +begin seed species + + # Start with receptors in different methylation states + ReceptorA(m~1,bind) RecA_count + ReceptorB(m~2,bind) RecB_count + +end seed species + +begin observables + Molecules RecA_m0 ReceptorA(m~0) + Molecules RecA_m1 ReceptorA(m~1) + Molecules RecA_m2 ReceptorA(m~2) + Molecules RecA_m3 ReceptorA(m~3) + + Molecules RecB_m0 ReceptorB(m~0) + Molecules RecB_m1 ReceptorB(m~1) + Molecules RecB_m2 ReceptorB(m~2) + Molecules RecB_m3 ReceptorB(m~3) + + Molecules Complex ReceptorA(bind!1).ReceptorB(bind!1) + +end observables + +begin functions + + # Separate local functions for each reactant type + # These will be combined in a DOR2 reaction + rateA(x) = k_base * (1 + RecA_m1(x) + 2*RecA_m2(x) + 3*RecA_m3(x)) + rateB(y) = k_base * (1 + RecB_m1(y) + 2*RecB_m2(y) + 3*RecB_m3(y)) + +end functions + +begin reaction rules + + # Bimolecular reaction with separate DOR functions for each reactant + # This should trigger DOR2RxnClass creation + ReceptorA%x(bind) + ReceptorB%y(bind) -> ReceptorA%x(bind!1).ReceptorB%y(bind!1) rateA(x)*rateB(y) + + # Unbinding reaction + ReceptorA(bind!1).ReceptorB(bind!1) -> ReceptorA(bind) + ReceptorB(bind) 0.1 + + # Some methylation/demethylation reactions to change states + ReceptorA(m~0) -> ReceptorA(m~1) 0.5 + ReceptorA(m~1) -> ReceptorA(m~2) 0.5 + ReceptorA(m~2) -> ReceptorA(m~3) 0.5 + ReceptorA(m~3) -> ReceptorA(m~2) 0.2 + ReceptorA(m~2) -> ReceptorA(m~1) 0.2 + ReceptorA(m~1) -> ReceptorA(m~0) 0.2 + + ReceptorB(m~0) -> ReceptorB(m~1) 0.5 + ReceptorB(m~1) -> ReceptorB(m~2) 0.5 + ReceptorB(m~2) -> ReceptorB(m~3) 0.5 + ReceptorB(m~3) -> ReceptorB(m~2) 0.2 + ReceptorB(m~2) -> ReceptorB(m~1) 0.2 + ReceptorB(m~1) -> ReceptorB(m~0) 0.2 + +end reaction rules + +writeXML(); \ No newline at end of file diff --git a/test/testSuite/t_dor2.xml b/test/testSuite/t_dor2.xml new file mode 100644 index 00000000..df931a3f --- /dev/null +++ b/test/testSuite/t_dor2.xml @@ -0,0 +1,789 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + k_base*(((1+RecA_m1(x))+(2*RecA_m2(x)))+(3*RecA_m3(x))) + + + + + + + + + + + + + k_base*(((1+RecB_m1(y))+(2*RecB_m2(y)))+(3*RecB_m3(y))) + + + + + + + + + + + + + rateA(x)*rateB(y) + + + + diff --git a/validate/basicModels/r18.txt b/validate/basicModels/r18.txt index 79e424db..e9f0161a 100644 --- a/validate/basicModels/r18.txt +++ b/validate/basicModels/r18.txt @@ -1,2 +1,2 @@ Testing for component handling with basic reactions in molecules with symmetric components --sim 5 -oSteps 100 \ No newline at end of file +-sim 10 -oSteps 100 -cb \ No newline at end of file diff --git a/validate/basicModels/r19.txt b/validate/basicModels/r19.txt new file mode 100644 index 00000000..d5ae7bc3 --- /dev/null +++ b/validate/basicModels/r19.txt @@ -0,0 +1,2 @@ +Testing for the correct handling of parameters labeled through species-scoped references +-sim 10 -oSteps 100 \ No newline at end of file diff --git a/validate/basicModels/r20.txt b/validate/basicModels/r20.txt new file mode 100644 index 00000000..94fb2d82 --- /dev/null +++ b/validate/basicModels/r20.txt @@ -0,0 +1,2 @@ +connectivity testing, simple binding, unbinding, and phosphorylation system +-sim 100 -oSteps 100 -connect -rxnlog basicModels/v20_rxns.json -trackconnected -printconnected -printmoltypes -printrxncounts diff --git a/validate/basicModels/r21.txt b/validate/basicModels/r21.txt new file mode 100644 index 00000000..f3338ca2 --- /dev/null +++ b/validate/basicModels/r21.txt @@ -0,0 +1,2 @@ +connectivity testing, mRNA monomer model +-sim 100 -oSteps 100 -connect -rxnlog basicModels/v21_rxns.json -utl 3 -gml 1000000 -trackconnected -printconnected -printmoltypes -printrxncounts diff --git a/validate/basicModels/r22.txt b/validate/basicModels/r22.txt new file mode 100644 index 00000000..1b08386c --- /dev/null +++ b/validate/basicModels/r22.txt @@ -0,0 +1,2 @@ +connectivity testing, mRNA monomer model with add operations +-sim 100 -oSteps 100 -connect -rxnlog basicModels/v22_rxns.json -utl 3 -gml 1000000 -trackconnected -printconnected -printmoltypes -printrxncounts diff --git a/validate/basicModels/r23.txt b/validate/basicModels/r23.txt new file mode 100644 index 00000000..90c34877 --- /dev/null +++ b/validate/basicModels/r23.txt @@ -0,0 +1,2 @@ +connectivity testing, mRNA monomer model with delete operations +-sim 100 -oSteps 100 -connect -rxnlog basicModels/v23_rxns.json -utl 3 -gml 1000000 -trackconnected -printconnected -printmoltypes -printrxncounts diff --git a/validate/basicModels/r24.txt b/validate/basicModels/r24.txt new file mode 100644 index 00000000..802a98c0 --- /dev/null +++ b/validate/basicModels/r24.txt @@ -0,0 +1,2 @@ +connectivity testing, mRNA monomer model with add and delete operations +-sim 100 -oSteps 100 -connect -rxnlog basicModels/v24_rxns.json -utl 3 -gml 1000000 -trackconnected -printconnected -printmoltypes -printrxncounts \ No newline at end of file diff --git a/validate/basicModels/r25.txt b/validate/basicModels/r25.txt new file mode 100644 index 00000000..c4c82e04 --- /dev/null +++ b/validate/basicModels/r25.txt @@ -0,0 +1,2 @@ +connectivity testing, mRNA monomer model, only `-rxnlog` check +-sim 100 -oSteps 100 -rxnlog basicModels/v25_rxns.json -utl 3 -gml 1000000 \ No newline at end of file diff --git a/validate/basicModels/r26.txt b/validate/basicModels/r26.txt new file mode 100644 index 00000000..814bf131 --- /dev/null +++ b/validate/basicModels/r26.txt @@ -0,0 +1,2 @@ +connectivity testing, mRNA monomer model, logbuffer size check +-sim 100 -oSteps 100 -rxnlog basicModels/v25_rxns.json -logbuffer 10000 -utl 3 -gml 1000000 \ No newline at end of file diff --git a/validate/basicModels/r27.txt b/validate/basicModels/r27.txt new file mode 100644 index 00000000..404782d9 --- /dev/null +++ b/validate/basicModels/r27.txt @@ -0,0 +1,2 @@ +NFSIM ONLY: population testing +-sim 100 -oSteps 100 \ No newline at end of file diff --git a/validate/basicModels/r28.txt b/validate/basicModels/r28.txt new file mode 100644 index 00000000..0c511e64 --- /dev/null +++ b/validate/basicModels/r28.txt @@ -0,0 +1,2 @@ +NFSIM ONLY: population testing, with reaction log +-sim 100 -oSteps 100 -rxnlog basicModels/v28_rxns.json \ No newline at end of file diff --git a/validate/basicModels/r29.txt b/validate/basicModels/r29.txt new file mode 100644 index 00000000..1cfa5b63 --- /dev/null +++ b/validate/basicModels/r29.txt @@ -0,0 +1,2 @@ +NFSIM ONLY: population testing, with reaction log +-sim 100 -oSteps 100 -rxnlog basicModels/v29_rxns.json \ No newline at end of file diff --git a/validate/basicModels/r30.txt b/validate/basicModels/r30.txt new file mode 100644 index 00000000..e18c552a --- /dev/null +++ b/validate/basicModels/r30.txt @@ -0,0 +1,2 @@ +NFSIM ONLY: population testing, with reaction log +-sim 100 -oSteps 100 -rxnlog basicModels/v30_rxns.json \ No newline at end of file diff --git a/validate/basicModels/v17.bngl b/validate/basicModels/v17.bngl index 7ac4242d..0ce817c8 100644 --- a/validate/basicModels/v17.bngl +++ b/validate/basicModels/v17.bngl @@ -14,7 +14,7 @@ R(s!1,s).L(s!2,s!1).R(s!2,s!3).L(s!3,s!4).R(s!4,s) 1 end seed species begin observables -Molecules Rpattern L(s!2,s!1).R(s!2,s!3).L(s!3,s!4) +Molecules Rpattern L(s!2,s!+).R(s!2,s!3).L(s!3,s!+) Molecules FreeR R(s,s) Molecules Rsimple R(s!1).L(s!1,s) Molecules Rtot R() diff --git a/validate/basicModels/v18.bngl b/validate/basicModels/v18.bngl index 47158688..0a6b8ee7 100644 --- a/validate/basicModels/v18.bngl +++ b/validate/basicModels/v18.bngl @@ -37,5 +37,6 @@ begin actions resetConcentrations() simulate({suffix=>"ssa",method=>"ssa",t_start=>0,t_end=>10,n_steps=>100}) resetConcentrations() - simulate({suffix=>"nf",method=>"nf",complex=>1,t_start=>0,t_end=>10,n_steps=>100}) -end actions \ No newline at end of file + writeXML() + #simulate({suffix=>"nf",method=>"nf",complex=>1,t_start=>0,t_end=>10,n_steps=>100}) +end actions diff --git a/validate/basicModels/v19.bngl b/validate/basicModels/v19.bngl new file mode 100644 index 00000000..7668753f --- /dev/null +++ b/validate/basicModels/v19.bngl @@ -0,0 +1,62 @@ +begin model + +begin parameters + + +kmod 1 + +end parameters + +begin molecule types + +X(a,a) +Y(b,c,d~0~1) + +A(a,a) +B(b,c,d~0~1) + +end molecule types + +begin seed species +X(a!1,a!2).Y(b!1,c,d~0).Y(b!2,c,d~0) 100 +A(a!1,a!2).B(b!1,c,d~0).B(b!2,c,d~0) 100 + +end seed species + +begin observables + +Molecules Ytot Y() +Molecules Yzero Y(d~0) + +Molecules Btot B() +Molecules Bzero B(d~0) + +end observables + +begin functions +factor(x) if(Ytot(x)>0,1,0) +factor2(x) if(Btot(x)>0,1,0) + +end functions + +begin reaction rules + +%z:X(a!1,a!+).Y(b!1,d~0)-> %z:X(a!1,a!+).Y(b!1,d~1) factor(z)*kmod # factor(z) should evaluate to 1 +%z:B(b!1,d~0).A(a!1,a!+)->%z:B(b!1,d~1).A(a!1,a!+) factor2(z)*kmod # factor(z) should evaluate to 1 + + +end reaction rules + +end model + + +generate_network({overwrite=>1}) +simulate({suffix=>"ode",method=>"ode",t_start=>0,t_end=>10,n_steps=>100}) +resetConcentrations() +simulate({suffix=>"ssa",method=>"ssa",t_start=>0,t_end=>10,n_steps=>100}) +resetConcentrations() +writeXML() + +#simulate({method=>"nf",gml=>2000000,complex=>1,get_final_state=>0,print_functions=>1,\ + t_start=>0,t_end=>30,n_steps=>300}) + diff --git a/validate/basicModels/v20.bngl b/validate/basicModels/v20.bngl new file mode 100644 index 00000000..e6e107b1 --- /dev/null +++ b/validate/basicModels/v20.bngl @@ -0,0 +1,32 @@ +# simple_system.bngl +begin model +begin parameters + kon 15 + koff 10 + kcat 0.6 + dephos 0.5 +end parameters +begin species + X(y,p~0) 5000 + X(y,p~1) 0 + Y(x) 500 +end species +begin reaction rules + X(p~1) -> X(p~0) dephos + X(y,p~0) + Y(x) -> X(y!1,p~0).Y(x!1) kon + X(y!1,p~0).Y(x!1) -> X(y,p~0) + Y(x) koff + X(y!1,p~0).Y(x!1) -> X(y,p~1) + Y(x) kcat +end reaction rules +begin observables + #Molecules X_free X(p~0,y) + Molecules Xp_free X(p~1,y) + #Molecules XY X(y!1).Y(x!1) +end observables +end model + +generate_network({overwrite=>1}) +simulate_ode({suffix=>"ode",t_end=>100,n_steps=>100}) +resetConcentrations() +simulate_ssa({suffix=>"ssa",t_end=>100,n_steps=>100}) +resetConcentrations() +writeXML() diff --git a/validate/basicModels/v21.bngl b/validate/basicModels/v21.bngl new file mode 100644 index 00000000..5fda3b53 --- /dev/null +++ b/validate/basicModels/v21.bngl @@ -0,0 +1,47 @@ +begin model +begin parameters + l_ribo 1 + l_mrna 10 + k_init 1 + k_elong 1 + k_term 1 + n_mrna_0 1 + n_ribo_0 100 +end parameters + +begin molecule types + mrna(p0~free~blocked,p1,p2,p3,p4,p5,p6,p7,p8,p9) + ribosome(asite) +end molecule types + +begin species + mrna(p0~free,p1,p2,p3,p4,p5,p6,p7,p8,p9) n_mrna_0 + ribosome(asite) n_ribo_0 +end species + +begin observables + Molecules rib_elong ribosome(asite!1).mrna(p1,p2!1) +end observables + +begin reaction rules + initiate: ribosome(asite) + mrna(p0~free) -> ribosome(asite!1).mrna(p0~blocked!1) k_init TotalRate + elongate_1: ribosome(asite!1).mrna(p0~blocked!1,p1) -> ribosome(asite!1).mrna(p0~free,p1!1) k_elong + elongate_2: ribosome(asite!1).mrna(p1!1,p2) -> ribosome(asite!1).mrna(p1,p2!1) k_elong + elongate_3: ribosome(asite!1).mrna(p2!1,p3) -> ribosome(asite!1).mrna(p2,p3!1) k_elong + elongate_4: ribosome(asite!1).mrna(p3!1,p4) -> ribosome(asite!1).mrna(p3,p4!1) k_elong + elongate_5: ribosome(asite!1).mrna(p4!1,p5) -> ribosome(asite!1).mrna(p4,p5!1) k_elong + elongate_6: ribosome(asite!1).mrna(p5!1,p6) -> ribosome(asite!1).mrna(p5,p6!1) k_elong + elongate_7: ribosome(asite!1).mrna(p6!1,p7) -> ribosome(asite!1).mrna(p6,p7!1) k_elong + elongate_8: ribosome(asite!1).mrna(p7!1,p8) -> ribosome(asite!1).mrna(p7,p8!1) k_elong + elongate_9: ribosome(asite!1).mrna(p8!1,p9) -> ribosome(asite!1).mrna(p8,p9!1) k_elong + terminate: ribosome(asite!1).mrna(p9!1) -> ribosome(asite) + mrna(p9) k_term +end reaction rules + +end model + +generate_network({overwrite=>1}) +simulate_ode({suffix=>"ode",t_end=>100,n_steps=>100}) +resetConcentrations() +simulate_ssa({suffix=>"ssa",t_end=>100,n_steps=>100}) +resetConcentrations() +writeXML() diff --git a/validate/basicModels/v22.bngl b/validate/basicModels/v22.bngl new file mode 100644 index 00000000..57b465d1 --- /dev/null +++ b/validate/basicModels/v22.bngl @@ -0,0 +1,45 @@ +begin model +begin parameters + l_ribo 1 + l_mrna 10 + add 0.01 + k_add 1 + k_init 1 + k_elong 1 + k_term 1 + n_mrna_0 1 + n_ribo_0 100 +end parameters + +begin molecule types + mrna(p0~free~blocked,p1,p2,p3) + ribosome(asite) +end molecule types + +begin species + mrna(p0~free,p1,p2,p3) n_mrna_0 + ribosome(asite) n_ribo_0 +end species + +begin observables + Molecules rib_elong ribosome(asite!1).mrna(p1,p2!1) +end observables + +begin reaction rules + test_add: 0 -> mrna(p0~free,p1,p2,p3) add + test_add_2: 0 -> ribosome(asite!1).mrna(p0~blocked!1,p1,p2,p3) add + initiate: ribosome(asite) + mrna(p0~free) -> ribosome(asite!1).mrna(p0~blocked!1) k_init TotalRate + elongate_1: ribosome(asite!1).mrna(p0~blocked!1,p1) -> ribosome(asite!1).mrna(p0~free,p1!1) k_elong + elongate_2: ribosome(asite!1).mrna(p1!1,p2) -> ribosome(asite!1).mrna(p1,p2!1) k_elong + elongate_3: ribosome(asite!1).mrna(p2!1,p3) -> ribosome(asite!1).mrna(p2,p3!1) k_elong + terminate: ribosome(asite!1).mrna(p3!1) -> ribosome(asite) + mrna(p3) k_term +end reaction rules + +end model + +generate_network({overwrite=>1}) +simulate_ode({suffix=>"ode",t_end=>100,n_steps=>100}) +resetConcentrations() +simulate_ssa({suffix=>"ssa",t_end=>100,n_steps=>100}) +resetConcentrations() +writeXML() \ No newline at end of file diff --git a/validate/basicModels/v23.bngl b/validate/basicModels/v23.bngl new file mode 100644 index 00000000..d8a86f4e --- /dev/null +++ b/validate/basicModels/v23.bngl @@ -0,0 +1,44 @@ +begin model +begin parameters + l_ribo 1 + l_mrna 10 + k_del 0.1 + k_init 1 + k_elong 1 + k_term 1 + n_mrna_0 1 + n_ribo_0 100 +end parameters + +begin molecule types + mrna(p0~free~blocked,p1,p2,p3) + ribosome(asite) +end molecule types + +begin species + mrna(p0~free,p1,p2,p3) n_mrna_0 + ribosome(asite) n_ribo_0 +end species + +begin observables + Molecules rib_elong ribosome(asite!1).mrna(p1,p2!1) +end observables + +begin reaction rules + test_del: ribosome(asite) -> 0 k_del + test_del_2: ribosome(asite!1).mrna(p0~blocked!1,p1,p2,p3) -> 0 k_del + initiate: ribosome(asite) + mrna(p0~free) -> ribosome(asite!1).mrna(p0~blocked!1) k_init TotalRate + elongate_1: ribosome(asite!1).mrna(p0~blocked!1,p1) -> ribosome(asite!1).mrna(p0~free,p1!1) k_elong + elongate_2: ribosome(asite!1).mrna(p1!1,p2) -> ribosome(asite!1).mrna(p1,p2!1) k_elong + elongate_3: ribosome(asite!1).mrna(p2!1,p3) -> ribosome(asite!1).mrna(p2,p3!1) k_elong + terminate: ribosome(asite!1).mrna(p3!1) -> ribosome(asite) + mrna(p3) k_term +end reaction rules + +end model + +generate_network({overwrite=>1}) +simulate_ode({suffix=>"ode",t_end=>100,n_steps=>100}) +resetConcentrations() +simulate_ssa({suffix=>"ssa",t_end=>100,n_steps=>100}) +resetConcentrations() +writeXML() diff --git a/validate/basicModels/v24.bngl b/validate/basicModels/v24.bngl new file mode 100644 index 00000000..98e85871 --- /dev/null +++ b/validate/basicModels/v24.bngl @@ -0,0 +1,47 @@ +begin model +begin parameters + l_ribo 1 + l_mrna 10 + k_add 1.2 + k_del 1 + k_init 1 + k_elong 1 + k_term 1 + n_mrna_0 1 + n_ribo_0 100 +end parameters + +begin molecule types + mrna(p0~free~blocked,p1,p2,p3) + ribosome(asite) +end molecule types + +begin species + mrna(p0~free,p1,p2,p3) n_mrna_0 + ribosome(asite) n_ribo_0 +end species + +begin observables + Molecules rib_elong ribosome(asite!1).mrna(p1,p2!1) +end observables + +begin reaction rules + test_add: 0 -> mrna(p0~free,p1,p2,p3) k_add + test_add_2: 0 -> ribosome(asite!1).mrna(p0~blocked!1,p1,p2,p3) k_add + test_del: ribosome(asite) -> 0 k_del + test_del_2: ribosome(asite!1).mrna(p0~blocked!1,p1,p2,p3) -> 0 k_del + initiate: ribosome(asite) + mrna(p0~free) -> ribosome(asite!1).mrna(p0~blocked!1) k_init TotalRate + elongate_1: ribosome(asite!1).mrna(p0~blocked!1,p1) -> ribosome(asite!1).mrna(p0~free,p1!1) k_elong + elongate_2: ribosome(asite!1).mrna(p1!1,p2) -> ribosome(asite!1).mrna(p1,p2!1) k_elong + elongate_3: ribosome(asite!1).mrna(p2!1,p3) -> ribosome(asite!1).mrna(p2,p3!1) k_elong + terminate: ribosome(asite!1).mrna(p3!1) -> ribosome(asite) + mrna(p3) k_term +end reaction rules + +end model + +generate_network({overwrite=>1}) +simulate_ode({suffix=>"ode",t_end=>100,n_steps=>100}) +resetConcentrations() +simulate_ssa({suffix=>"ssa",t_end=>100,n_steps=>100}) +resetConcentrations() +writeXML() diff --git a/validate/basicModels/v25.bngl b/validate/basicModels/v25.bngl new file mode 100644 index 00000000..98e85871 --- /dev/null +++ b/validate/basicModels/v25.bngl @@ -0,0 +1,47 @@ +begin model +begin parameters + l_ribo 1 + l_mrna 10 + k_add 1.2 + k_del 1 + k_init 1 + k_elong 1 + k_term 1 + n_mrna_0 1 + n_ribo_0 100 +end parameters + +begin molecule types + mrna(p0~free~blocked,p1,p2,p3) + ribosome(asite) +end molecule types + +begin species + mrna(p0~free,p1,p2,p3) n_mrna_0 + ribosome(asite) n_ribo_0 +end species + +begin observables + Molecules rib_elong ribosome(asite!1).mrna(p1,p2!1) +end observables + +begin reaction rules + test_add: 0 -> mrna(p0~free,p1,p2,p3) k_add + test_add_2: 0 -> ribosome(asite!1).mrna(p0~blocked!1,p1,p2,p3) k_add + test_del: ribosome(asite) -> 0 k_del + test_del_2: ribosome(asite!1).mrna(p0~blocked!1,p1,p2,p3) -> 0 k_del + initiate: ribosome(asite) + mrna(p0~free) -> ribosome(asite!1).mrna(p0~blocked!1) k_init TotalRate + elongate_1: ribosome(asite!1).mrna(p0~blocked!1,p1) -> ribosome(asite!1).mrna(p0~free,p1!1) k_elong + elongate_2: ribosome(asite!1).mrna(p1!1,p2) -> ribosome(asite!1).mrna(p1,p2!1) k_elong + elongate_3: ribosome(asite!1).mrna(p2!1,p3) -> ribosome(asite!1).mrna(p2,p3!1) k_elong + terminate: ribosome(asite!1).mrna(p3!1) -> ribosome(asite) + mrna(p3) k_term +end reaction rules + +end model + +generate_network({overwrite=>1}) +simulate_ode({suffix=>"ode",t_end=>100,n_steps=>100}) +resetConcentrations() +simulate_ssa({suffix=>"ssa",t_end=>100,n_steps=>100}) +resetConcentrations() +writeXML() diff --git a/validate/basicModels/v26.bngl b/validate/basicModels/v26.bngl new file mode 100644 index 00000000..98e85871 --- /dev/null +++ b/validate/basicModels/v26.bngl @@ -0,0 +1,47 @@ +begin model +begin parameters + l_ribo 1 + l_mrna 10 + k_add 1.2 + k_del 1 + k_init 1 + k_elong 1 + k_term 1 + n_mrna_0 1 + n_ribo_0 100 +end parameters + +begin molecule types + mrna(p0~free~blocked,p1,p2,p3) + ribosome(asite) +end molecule types + +begin species + mrna(p0~free,p1,p2,p3) n_mrna_0 + ribosome(asite) n_ribo_0 +end species + +begin observables + Molecules rib_elong ribosome(asite!1).mrna(p1,p2!1) +end observables + +begin reaction rules + test_add: 0 -> mrna(p0~free,p1,p2,p3) k_add + test_add_2: 0 -> ribosome(asite!1).mrna(p0~blocked!1,p1,p2,p3) k_add + test_del: ribosome(asite) -> 0 k_del + test_del_2: ribosome(asite!1).mrna(p0~blocked!1,p1,p2,p3) -> 0 k_del + initiate: ribosome(asite) + mrna(p0~free) -> ribosome(asite!1).mrna(p0~blocked!1) k_init TotalRate + elongate_1: ribosome(asite!1).mrna(p0~blocked!1,p1) -> ribosome(asite!1).mrna(p0~free,p1!1) k_elong + elongate_2: ribosome(asite!1).mrna(p1!1,p2) -> ribosome(asite!1).mrna(p1,p2!1) k_elong + elongate_3: ribosome(asite!1).mrna(p2!1,p3) -> ribosome(asite!1).mrna(p2,p3!1) k_elong + terminate: ribosome(asite!1).mrna(p3!1) -> ribosome(asite) + mrna(p3) k_term +end reaction rules + +end model + +generate_network({overwrite=>1}) +simulate_ode({suffix=>"ode",t_end=>100,n_steps=>100}) +resetConcentrations() +simulate_ssa({suffix=>"ssa",t_end=>100,n_steps=>100}) +resetConcentrations() +writeXML() diff --git a/validate/basicModels/v27.bngl b/validate/basicModels/v27.bngl new file mode 100644 index 00000000..9fc89656 --- /dev/null +++ b/validate/basicModels/v27.bngl @@ -0,0 +1,114 @@ +# Created by BioNetGen 2.8.4 + +begin model +begin parameters + f 0.01 + NA 6.0221e23 + V 1e-12*f + L0 (500e-9*NA)*V + R0 (100e-9*NA)*V + A0 (100e-9*NA)*V + B0 (100e-9*NA)*V + C0 (100e-9*NA)*V + BC0 (100e-9*NA)*V + kp1 10e6/(NA*V) + km1 1.0 + k2 1.0 + k3 1.0 + kp4 10e6/(NA*V) + km4 1.0 + kp5 10e6/(NA*V) + km5 1.0 + kp6 10e6/(NA*V) + km6 1.0 + k7 1.0 + k8 1.0 + kp9 10e6/(NA*V) + km9 1.0 + k_lump 10000 +end parameters +begin molecule types + A(r,b~0~P) + B(r,c) + C(b) + L(r) + P1() population + P2() population + P3() population + P4() population + P5() population + P6() population + P7() population + P8() population + R(l,a~0~P,b~0~P) +end molecule types +begin observables + Molecules LR L(r!1).R(l!1) + Molecules Rp R(a~P!?) R(b~P!?) + Molecules Ap A(b~P!?) P3() P4() P5() + Molecules RC R(b~P!1).B(r!1,c!2).C(b!2) R(a~P!1).A(r!1,b~P!2).B(r!2,c!3).C(b!3) + Molecules Ltot L() P1() + Molecules Rtot R() + Molecules Atot A() P2() P3() P4() P5() + Molecules Btot B() P4() P5() P6() P7() + Molecules Ctot C() P5() P7() P8() +end observables +begin species + P1() L0 + R(a~0,b~0,l) R0 + P2() A0 + P6() B0 + P8() C0 + P7() BC0 + P3() 0 + P4() 0 + P5() 0 +end species +begin reaction rules + _R1_v1: P1() + R(l) -> L(r!1).R(l!1) kp1 + _reverse__R1_v1: L(r!1).R(l!1) -> P1() + R(l) km1 DeleteMolecules + _R2_v1: L(r!1).R(a~0,l!1) -> L(r!1).R(a~P,l!1) k2 + _R3_v1: L(r!1).R(b~0,l!1) -> L(r!1).R(b~P,l!1) k2 + _R4_v1: R(a~P) -> R(a~0) k3 + _R5_v1: R(b~P) -> R(b~0) k3 + _R6_v1: R(a~P) + A(r) -> A(r!1).R(a~P!1) kp4 + _R6_v2: R(a~P) + P2() -> A(b~0,r!1).R(a~P!1) kp4 + _R6_v3: R(a~P) + P3() -> A(b~P,r!1).R(a~P!1) kp4 + _R6_v4: R(a~P) + P4() -> A(b~P!1,r!2).B(c,r!1).R(a~P!2) kp4 + _R6_v5: R(a~P) + P5() -> A(b~P!1,r!2).B(c!3,r!1).C(b!3).R(a~P!2) kp4 + _reverse__R6_v1: A(r!1).R(a~P!1) -> R(a~P) + A(r) km4 + _R7_v1: R(b~P) + B(r) -> B(r!1).R(b~P!1) kp5 + _R7_v2: R(b~P) + P6() -> B(c,r!1).R(b~P!1) kp5 + _R7_v3: R(b~P) + P7() -> B(c!1,r!2).C(b!1).R(b~P!2) kp5 + _reverse__R7_v1: B(r!1).R(b~P!1) -> R(b~P) + B(r) km5 + _R8_v1: B(c) + P8() -> B(c!1).C(b!1) kp6 + _R8_v2: P4() + P8() -> P5() kp6 + _R8_v3: P6() + P8() -> P7() kp6 + _reverse__R8_v1: B(c!1).C(b!1) -> B(c) + P8() km6 DeleteMolecules + _reverse__R8_v2: P5() -> P4() + P8() km6 + _reverse__R8_v3: P7() -> P6() + P8() km6 + _R9_v1: A(b~0,r!1).R(a~P!1) -> A(b~P,r!1).R(a~P!1) k7 + _R10_v1: A(b~P) -> A(b~0) k8 + _R10_v2: P3() -> P2() k8 + _R11_v1: A(b~P) + B(r) -> A(b~P!1).B(r!1) kp9 + _R11_v2: A(b~P) + P6() -> A(b~P!1).B(c,r!1) kp9 + _R11_v3: A(b~P) + P7() -> A(b~P!1).B(c!2,r!1).C(b!2) kp9 + _R11_v4: P3() + B(r) -> A(b~P!1,r).B(r!1) kp9 + _R11_v5: P3() + P6() -> P4() kp9 + _R11_v6: P3() + P7() -> P5() kp9 + _reverse__R11_v1: A(b~P!1).B(r!1) -> A(b~P) + B(r) km9 + _reverse__R11_v2: P4() -> P3() + P6() km9 + _reverse__R11_v3: P5() -> P3() + P7() km9 + _MapRule0: L(r) -> P1() k_lump + _MapRule1: A(b~0,r) -> P2() k_lump + _MapRule2: A(b~P,r) -> P3() k_lump + _MapRule3: A(b~P!1,r).B(c,r!1) -> P4() k_lump + _MapRule4: A(b~P!1,r).B(c!2,r!1).C(b!2) -> P5() k_lump + _MapRule5: B(c,r) -> P6() k_lump + _MapRule6: B(c!1,r).C(b!1) -> P7() k_lump + _MapRule7: C(b) -> P8() k_lump +end reaction rules +end model + +## model actions ## +writeXML() diff --git a/validate/basicModels/v28.bngl b/validate/basicModels/v28.bngl new file mode 100644 index 00000000..47c1ca4b --- /dev/null +++ b/validate/basicModels/v28.bngl @@ -0,0 +1,116 @@ +# Created by BioNetGen 2.8.4 + +begin model +begin parameters + f 0.01 + NA 6.0221e23 + V 1e-12*f + L0 (500e-9*NA)*V + R0 (100e-9*NA)*V + A0 (100e-9*NA)*V + B0 (100e-9*NA)*V + C0 (100e-9*NA)*V + BC0 (100e-9*NA)*V + kp1 10e6/(NA*V) + km1 1.0 + k2 1.0 + k3 1.0 + kp4 10e6/(NA*V) + km4 1.0 + kp5 10e6/(NA*V) + km5 1.0 + kp6 10e6/(NA*V) + km6 1.0 + k7 1.0 + k8 1.0 + kp9 10e6/(NA*V) + km9 1.0 + k_lump 10000 +end parameters +begin molecule types + A(r,b~0~P) + B(r,c) + C(b) + L(r) + P1() population + P2() population + P3() population + P4() population + P5() population + P6() population + P7() population + P8() population + R(l,a~0~P,b~0~P) +end molecule types +begin observables + Molecules LR L(r!1).R(l!1) + Molecules Rp R(a~P!?) R(b~P!?) + Molecules Ap A(b~P!?) P3() P4() P5() + Molecules RC R(b~P!1).B(r!1,c!2).C(b!2) R(a~P!1).A(r!1,b~P!2).B(r!2,c!3).C(b!3) + Molecules Ltot L() P1() + Molecules Rtot R() + Molecules Atot A() P2() P3() P4() P5() + Molecules Btot B() P4() P5() P6() P7() + Molecules Ctot C() P5() P7() P8() +end observables +begin species + P1() L0 + R(a~0,b~0,l) R0 + P2() A0 + P6() B0 + P8() C0 + P7() BC0 + P3() 0 + P4() 0 + P5() 0 +end species +begin reaction rules + _R1_v1: P1() + R(l) -> L(r!1).R(l!1) kp1 + _reverse__R1_v1: L(r!1).R(l!1) -> P1() + R(l) km1 DeleteMolecules + _R2_v1: L(r!1).R(a~0,l!1) -> L(r!1).R(a~P,l!1) k2 + _R3_v1: L(r!1).R(b~0,l!1) -> L(r!1).R(b~P,l!1) k2 + _R4_v1: R(a~P) -> R(a~0) k3 + _R5_v1: R(b~P) -> R(b~0) k3 + _R6_v1: R(a~P) + A(r) -> A(r!1).R(a~P!1) kp4 + _R6_v2: R(a~P) + P2() -> A(b~0,r!1).R(a~P!1) kp4 + _R6_v3: R(a~P) + P3() -> A(b~P,r!1).R(a~P!1) kp4 + _R6_v4: R(a~P) + P4() -> A(b~P!1,r!2).B(c,r!1).R(a~P!2) kp4 + _R6_v5: R(a~P) + P5() -> A(b~P!1,r!2).B(c!3,r!1).C(b!3).R(a~P!2) kp4 + _reverse__R6_v1: A(r!1).R(a~P!1) -> R(a~P) + A(r) km4 + _R7_v1: R(b~P) + B(r) -> B(r!1).R(b~P!1) kp5 + _R7_v2: R(b~P) + P6() -> B(c,r!1).R(b~P!1) kp5 + _R7_v3: R(b~P) + P7() -> B(c!1,r!2).C(b!1).R(b~P!2) kp5 + _reverse__R7_v1: B(r!1).R(b~P!1) -> R(b~P) + B(r) km5 + _R8_v1: B(c) + P8() -> B(c!1).C(b!1) kp6 + _R8_v2: P4() + P8() -> P5() kp6 + _R8_v3: P6() + P8() -> P7() kp6 + _reverse__R8_v1: B(c!1).C(b!1) -> B(c) + P8() km6 DeleteMolecules + _reverse__R8_v2: P5() -> P4() + P8() km6 + _reverse__R8_v3: P7() -> P6() + P8() km6 + _R9_v1: A(b~0,r!1).R(a~P!1) -> A(b~P,r!1).R(a~P!1) k7 + _R10_v1: A(b~P) -> A(b~0) k8 + _R10_v2: P3() -> P2() k8 + _R11_v1: A(b~P) + B(r) -> A(b~P!1).B(r!1) kp9 + _R11_v2: A(b~P) + P6() -> A(b~P!1).B(c,r!1) kp9 + _R11_v3: A(b~P) + P7() -> A(b~P!1).B(c!2,r!1).C(b!2) kp9 + _R11_v4: P3() + B(r) -> A(b~P!1,r).B(r!1) kp9 + _R11_v5: P3() + P6() -> P4() kp9 + _R11_v6: P3() + P7() -> P5() kp9 + _reverse__R11_v1: A(b~P!1).B(r!1) -> A(b~P) + B(r) km9 + _reverse__R11_v2: P4() -> P3() + P6() km9 + _reverse__R11_v3: P5() -> P3() + P7() km9 + _MapRule0: L(r) -> P1() k_lump + _MapRule1: A(b~0,r) -> P2() k_lump + _MapRule2: A(b~P,r) -> P3() k_lump + _MapRule3: A(b~P!1,r).B(c,r!1) -> P4() k_lump + _MapRule4: A(b~P!1,r).B(c!2,r!1).C(b!2) -> P5() k_lump + _MapRule5: B(c,r) -> P6() k_lump + _MapRule6: B(c!1,r).C(b!1) -> P7() k_lump + _MapRule7: C(b) -> P8() k_lump +end reaction rules +end model + + +## model actions ## +writeXML() + diff --git a/validate/basicModels/v29.bngl b/validate/basicModels/v29.bngl new file mode 100644 index 00000000..be07a2a9 --- /dev/null +++ b/validate/basicModels/v29.bngl @@ -0,0 +1,16 @@ +# issue37_test.bngl + +begin model +begin species + A() 0 +end species +begin observables + Molecules A A() +end observables +begin reaction rules + 0 <-> A() 10,1 +end reaction rules +end model + +## model actions ## +writeXML() \ No newline at end of file diff --git a/validate/basicModels/v30.bngl b/validate/basicModels/v30.bngl new file mode 100644 index 00000000..dcfdadf2 --- /dev/null +++ b/validate/basicModels/v30.bngl @@ -0,0 +1,17 @@ +# issue37_test2.bngl + +begin model +begin species + A() 100 + B() 0 +end species +begin observables + Molecules A A() +end observables +begin reaction rules + 0 <-> A() 10,1 +end reaction rules +end model + +## model actions ## +writeXML() \ No newline at end of file diff --git a/validate/requirements.txt b/validate/requirements.txt new file mode 100644 index 00000000..3d646db7 --- /dev/null +++ b/validate/requirements.txt @@ -0,0 +1,3 @@ +numpy +bionetgen +seaborn \ No newline at end of file diff --git a/validate/validate.py b/validate/validate.py new file mode 100644 index 00000000..2bf63c7e --- /dev/null +++ b/validate/validate.py @@ -0,0 +1,144 @@ +import unittest +import os +import numpy as np +import subprocess +import re +import fnmatch +import sys +import bionetgen + +nIterations=15 +nfsimPrePath='..' +mfolder='./basicModels' +bngPath = os.path.join(bionetgen.defaults.bng_path, "BNG2.pl") +if os.name == "nt": + nfsimPath = os.path.join(nfsimPrePath, 'build', 'NFsim.exe') +else: + nfsimPath = os.path.join(nfsimPrePath, 'build', 'NFsim') + + + +class ParametrizedTestCase(unittest.TestCase): + + """ TestCase classes that want to be parametrized should + inherit from this class. + """ + + def __init__(self, methodName='runTest', param=None): + super(ParametrizedTestCase, self).__init__(methodName) + self.param = param + + @staticmethod + def parametrize(testcase_klass, param=None): + """ Create a suite containing all tests taken from the given + subclass, passing them the parameter 'param'. + """ + testloader = unittest.TestLoader() + testnames = testloader.getTestCaseNames(testcase_klass) + suite = unittest.TestSuite() + for name in testnames: + suite.addTest(testcase_klass(name, param=param)) + return suite + + +def loadResults(fileName, split): + try: + with open(fileName) as dataInput: + timeCourse = [] + # remove spaces + line = dataInput.readline().strip() + headers = re.sub('\s+', ' ', line).split(split) + + for line in dataInput: + nline = re.sub('\s+', ' ', line.strip()).split(' ') + try: + timeCourse.append([float(x) for x in nline]) + except: + print('++++', nline) + return headers, np.array(timeCourse) + except IOError: + print('no file') + return [], [] + + +class TestNFSimFile(ParametrizedTestCase): + # XXX:ideally this should be done through the console but I'm doing the quick and dirty version right now + def BNGtrajectoryGeneration(self, outputDirectory, fileNumber): + bngFileName = os.path.join(outputDirectory, 'v{0}.bngl'.format(fileNumber)) + with open(os.devnull, "w") as fnull: + subprocess.check_call(['perl', bngPath, '-outdir', outputDirectory, '-log', bngFileName], stdout=fnull) + + def NFsimtrajectoryGeneration(self, outputDirectory, fileNumber, runOptions): + runOptions = [x.strip() for x in runOptions.split(' ')] + with open(os.devnull, "w") as fnull: + subprocess.check_call([nfsimPath, '-xml', os.path.join(outputDirectory, 'v{0}.xml'.format(fileNumber)), + '-o', os.path.join(outputDirectory, 'v{0}_nf.gdat'.format(fileNumber))] + runOptions, + stdout=fnull) + + def loadConfigurationFile(self, outputDirectory, fileNumber): + with open(os.path.join(outputDirectory, 'r{0}.txt').format(fileNumber), 'r') as f: + return f.readlines() + + def test_nfsim(self): + tol = 0.35 # this is the error tolerance when comparing nfsim's run to the ssa where 0.35 = 35% + (modelName, runOptions) = self.loadConfigurationFile(self.param['odir'], self.param['num']) + print(f"Processing model r{self.param['num']}.txt: {modelName.strip()}") + # here we decide if this is a NFsim only run or not + if modelName.startswith("NFSIM ONLY"): + self.BNGtrajectoryGeneration(self.param['odir'], self.param['num']) + self.NFsimtrajectoryGeneration(self.param['odir'], self.param['num'], runOptions) + nfh, nf = loadResults(os.path.join(self.param['odir'], 'v{0}_nf.gdat'.format(self.param['num'])), ' ') + # here we just need to make sure we managed to get here without errors + #assert len(nf) > 0 + self.assertTrue(len(nf) > 0) + else: + ssaDiff = nfDiff = 0 + for index in range(self.param['iterations']): + print(f'Iteration {index+1}') + self.BNGtrajectoryGeneration(self.param['odir'], self.param['num']) + self.NFsimtrajectoryGeneration(self.param['odir'], self.param['num'], runOptions) + odeh, ode = loadResults(os.path.join(self.param['odir'], 'v{0}_ode.gdat'.format(self.param['num'])), ' ') + ssah, ssa = loadResults(os.path.join(self.param['odir'], 'v{0}_ssa.gdat'.format(self.param['num'])), ' ') + nfh, nf = loadResults(os.path.join(self.param['odir'], 'v{0}_nf.gdat'.format(self.param['num'])), ' ') + + #square root difference + ssaDiff += pow(sum(pow(ode[:, 1:] - ssa[:, 1:], 2)), 0.5) + nfDiff += pow(sum(pow(ode[:, 1:] - nf[:, 1:], 2)), 0.5) + rdiff = nfDiff - ssaDiff - (tol * ssaDiff) + # relative difference should be less than 'tol' + bad=np.where(rdiff>0)[0] + if (bad.size>0): + print(f"Sir, the observables {bad+1} did not pass. Trying again") + else: + print("Check passed.") + break + self.assertTrue(bad.size==0) + + +def getTests(directory): + """ + Gets a list of bngl files that could be correctly translated in a given 'directory' + """ + matches = [] + for root, dirnames, filenames in os.walk(directory): + for filename in fnmatch.filter(filenames, '*txt'): + matches.append(''.join(filename.split('.')[0][1:])) + return sorted(matches) + +if __name__ == "__main__": + suite = unittest.TestSuite() + if len(sys.argv) > 1: + os.chdir(sys.argv[1]) + testFolder = mfolder + tests = getTests(testFolder) + for index in tests: + suite.addTest(ParametrizedTestCase.parametrize(TestNFSimFile, param={'num': index, + 'odir': mfolder, 'iterations': nIterations})) + result = unittest.TextTestRunner(verbosity=1).run(suite) + + ret = (list(result.failures) == [] and list(result.errors) == []) + ret = 0 if ret else 1 + if ret > 0: + sys.exit("Validation return an error code") + else: + sys.exit()