diff --git a/CMakeLists.txt b/CMakeLists.txt index 26d34a0444..107bb71708 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -8,7 +8,7 @@ cmake_minimum_required(VERSION 3.16) # Version information and include modules set(ROADRUNNER_VERSION_MAJOR 2) -set(ROADRUNNER_VERSION_MINOR 5) +set(ROADRUNNER_VERSION_MINOR 6) set(ROADRUNNER_VERSION_PATCH 0) set(ROADRUNNER_VERSION "${ROADRUNNER_VERSION_MAJOR}.${ROADRUNNER_VERSION_MINOR}.${ROADRUNNER_VERSION_PATCH}") diff --git a/azure-pipelines.yml b/azure-pipelines.yml index 942e933848..29eb8a3244 100644 --- a/azure-pipelines.yml +++ b/azure-pipelines.yml @@ -32,16 +32,30 @@ stages: timeoutInMinutes: "0" displayName: MacBuildRoadrunnerCpp continueOnError: "false" - pool: - vmImage: 'macOS-11' strategy: matrix: - 64-bit Mac Release: + 64-bit Mac Release x86: BuildType: Release LLVM_DOWNLOAD_LINK: 'https://github.com/sys-bio/llvm-13.x/releases/download/llvmorg-13.0.0/llvm-13.x-macosx_11_7_x86_64.zip' - 64-bit Mac Debug: + vmImage: macOS-11 + arch: x86_64 + 64-bit Mac Debug x86: BuildType: Debug LLVM_DOWNLOAD_LINK: 'https://github.com/sys-bio/llvm-13.x/releases/download/llvmorg-13.0.0/llvm-13.x-macosx_11_7_x86_64.zip' + vmImage: macOS-11 + arch: x86_64 + # 64-bit Mac Release M1: + # BuildType: Release + # LLVM_DOWNLOAD_LINK: 'https://github.com/sys-bio/llvm-13.x/releases/download/llvmorg-13.0.0/llvm-13.x-macosx_13_2_arm64.zip' + # vmImage: macOS-14 + # arch: arm64 + # 64-bit Mac Debug M1: + # BuildType: Debug + # LLVM_DOWNLOAD_LINK: 'https://github.com/sys-bio/llvm-13.x/releases/download/llvmorg-13.0.0/llvm-13.x-macosx_13_2_arm64.zip' + # vmImage: macOS-14 + # arch: arm64 + pool: + vmImage: $(vmImage) variables: LLVM_CACHE: 'false' PythonName: 'py39' @@ -138,7 +152,7 @@ stages: # - https://dev.azure.com/TheRoadrunnerProject/libroadrunner-deps/_build?definitionId=9 pipeline: 9 runVersion: 'latest' - artifact: libroadrunner-deps-$(Agent.OS)-$(BuildType) + artifact: libroadrunner-deps-$(Agent.OS)-$(BuildType)-$(arch) path: $(DEPS_INSTALL_PREFIX) displayName: Download libroadrunner-deps install artifacts - script: | @@ -191,22 +205,58 @@ stages: timeoutInMinutes: "0" displayName: MacBuildRoadrunnerPython continueOnError: "false" - pool: - vmImage: 'macOS-11' strategy: matrix: - py38: - PythonVersion: 3.8 - PythonName: py38 - py39: + py39, x86: + LLVM_DOWNLOAD_LINK: 'https://github.com/sys-bio/llvm-13.x/releases/download/llvmorg-13.0.0/llvm-13.x-macosx_11_7_x86_64.zip' + vmImage: macOS-11 + arch: x86_64 PythonVersion: 3.9 PythonName: py39 - py310: + py310, x86: + LLVM_DOWNLOAD_LINK: 'https://github.com/sys-bio/llvm-13.x/releases/download/llvmorg-13.0.0/llvm-13.x-macosx_11_7_x86_64.zip' + vmImage: macOS-11 + arch: x86_64 PythonVersion: 3.10 PythonName: py310 - py311: + py311, x86: + LLVM_DOWNLOAD_LINK: 'https://github.com/sys-bio/llvm-13.x/releases/download/llvmorg-13.0.0/llvm-13.x-macosx_11_7_x86_64.zip' + vmImage: macOS-11 + arch: x86_64 PythonVersion: 3.11 PythonName: py311 + py312, x86: + LLVM_DOWNLOAD_LINK: 'https://github.com/sys-bio/llvm-13.x/releases/download/llvmorg-13.0.0/llvm-13.x-macosx_11_7_x86_64.zip' + vmImage: macOS-11 + arch: x86_64 + PythonVersion: 3.12 + PythonName: py312 + # py39, arm64: + # LLVM_DOWNLOAD_LINK: 'https://github.com/sys-bio/llvm-13.x/releases/download/llvmorg-13.0.0/llvm-13.x-macosx_13_2_arm64.zip' + # vmImage: macOS-11 + # arch: arm64_64 + # PythonVersion: 3.9 + # PythonName: py39 + # py310, arm64: + # LLVM_DOWNLOAD_LINK: 'https://github.com/sys-bio/llvm-13.x/releases/download/llvmorg-13.0.0/llvm-13.x-macosx_13_2_arm64.zip' + # vmImage: macOS-11 + # arch: arm64_64 + # PythonVersion: 3.10 + # PythonName: py310 + # py311, arm64: + # LLVM_DOWNLOAD_LINK: 'https://github.com/sys-bio/llvm-13.x/releases/download/llvmorg-13.0.0/llvm-13.x-macosx_13_2_arm64.zip' + # vmImage: macOS-11 + # arch: arm64_64 + # PythonVersion: 3.11 + # PythonName: py311 + # py312, arm64: + # LLVM_DOWNLOAD_LINK: 'https://github.com/sys-bio/llvm-13.x/releases/download/llvmorg-13.0.0/llvm-13.x-macosx_13_2_arm64.zip' + # vmImage: macOS-11 + # arch: arm64_64 + # PythonVersion: 3.12 + # PythonName: py312 + pool: + vmImage: $(vmImage) variables: MinicondaRoot : '/usr/local/miniconda' PythonRoot: '$(MinicondaRoot)/envs/$(PythonName)' @@ -217,7 +267,6 @@ stages: SWIG_CACHE: 'false' MINICONDA_CACHE: 'false' LLVM_CACHE: 'false' - LLVM_DOWNLOAD_LINK: 'https://github.com/sys-bio/llvm-13.x/releases/download/llvmorg-13.0.0/llvm-13.x-macosx_11_7_x86_64.zip' steps: - checkout: self submodules: recursive @@ -272,7 +321,7 @@ stages: project: 'libroadrunner-deps' pipeline: 9 runVersion: 'latest' - artifact: libroadrunner-deps-$(Agent.OS)-Release + artifact: libroadrunner-deps-$(Agent.OS)-Release-$(arch) path: $(DEPS_INSTALL_PREFIX) displayName: Download libroadrunner-deps install artifacts - script: | @@ -370,6 +419,7 @@ stages: ls echo "$(PythonExecutable) setup.py bdist_wheel" $(PythonExecutable) setup.py bdist_wheel + rm -r build/ echo "$(PythonExecutable) setup_rrplugins.py bdist_wheel" $(PythonExecutable) setup_rrplugins.py bdist_wheel displayName: Generate pip wheel @@ -413,8 +463,8 @@ stages: variables: LLVM_CACHE: 'false' MinicondaRoot : 'C:\Miniconda' - PythonName: 'py39' - PythonVersion: '3.9' + PythonName: 'py311' + PythonVersion: '3.11' PythonRoot: '$(MinicondaRoot)\envs\$(PythonName)' PythonLibDir: '$(PythonRoot)\Lib' PythonScriptsDir: '$(PythonRoot)\Scripts' @@ -557,9 +607,6 @@ stages: vmImage: 'windows-2019' strategy: matrix: - py38: - PythonVersion: 3.8 - PythonName: py38 py39: PythonVersion: 3.9 PythonName: py39 @@ -569,6 +616,9 @@ stages: py311: PythonVersion: 3.11 PythonName: py311 + py312: + PythonVersion: 3.12 + PythonName: py312 variables: MinicondaRoot : 'C:\Miniconda' PythonRoot: '$(MinicondaRoot)\envs\$(PythonName)' @@ -733,6 +783,7 @@ stages: echo "ls in install dir" ls $(PythonExecutable) setup.py bdist_wheel + rm -r build/ $(PythonExecutable) setup_rrplugins.py bdist_wheel displayName: Generate pip wheel @@ -761,7 +812,7 @@ stages: displayName: UbuntuBuildRoadrunnerCpp continueOnError: "false" pool: - vmImage: 'Ubuntu-20.04' + vmImage: 'Ubuntu-latest' strategy: matrix: 64-bit Linux Release: @@ -963,9 +1014,6 @@ stages: - job: strategy: matrix: - py38: - PythonVersion: 3.8 - PythonName: py38 py39: PythonVersion: 3.9 PythonName: py39 @@ -975,8 +1023,11 @@ stages: py311: PythonVersion: 3.11 PythonName: py311 + py312: + PythonVersion: 3.12 + PythonName: py312 pool: - vmImage: 'ubuntu-20.04' + vmImage: 'ubuntu-latest' container: sysbiouw/roadrunner-manylinux2014:latest variables: CCACHE_DIR: '$(Pipeline.Workspace)/ccache' @@ -1074,6 +1125,7 @@ stages: echo "cmake command: cmake -DLLVM_INSTALL_PREFIX=/install-llvm-13.x -DRR_DEPENDENCIES_INSTALL_PREFIX=/install-libroadrunner-deps -DCMAKE_INSTALL_PREFIX=$(INSTALL_DIRECTORY) -DBUILD_PYTHON=ON -DBUILD_RR_PLUGINS=ON -DBUILD_TESTS=ON -DPython_ROOT_DIR=$(PythonRoot) -DSWIG_EXECUTABLE=$(SwigExecutable) -DCMAKE_BUILD_TYPE=Release .." cmake -DLLVM_INSTALL_PREFIX=/install-llvm-13.x -DRR_DEPENDENCIES_INSTALL_PREFIX=/install-libroadrunner-deps -DCMAKE_INSTALL_PREFIX=$(INSTALL_DIRECTORY) -DBUILD_PYTHON=ON -DBUILD_RR_PLUGINS=ON -DBUILD_TESTS=ON -DPython_ROOT_DIR=$(PythonRoot) -DSWIG_EXECUTABLE=$(SwigExecutable) -DCMAKE_BUILD_TYPE=Release .. + rm -r -f $(SOURCE_DIR)/docs $(SOURCE_DIR)/examples cmake --build . --target install --config Release -j 12 displayName: Build With Python @@ -1086,26 +1138,34 @@ stages: cd $(BUILD_DIR) ctest --extra-verbose --output-on-failure --exclude-regex python_tests_RunStochasticTestSuite . displayName: Run Tests - - script: | - echo ". ~/.bashrc" - . ~/.bashrc - $(CondaExecutable) activate $(PythonName) + # Skip stochastic test suite for now because it keeps pushing us over an hour. + # - script: | + # echo ". ~/.bashrc" + # . ~/.bashrc + # $(CondaExecutable) activate $(PythonName) - echo "cd'ing to build dir $(BUILD_DIR)" - cd $(BUILD_DIR) - ctest --extra-verbose --output-on-failure --tests-regex python_tests_RunStochasticTestSuite . - displayName: Run Stochastic test suite + # echo "cd'ing to build dir $(BUILD_DIR)" + # cd $(BUILD_DIR) + # ctest --extra-verbose --output-on-failure --tests-regex python_tests_RunStochasticTestSuite . + # displayName: Run Stochastic test suite - script: | echo "cd to $(INSTALL_DIRECTORY)" cd $(INSTALL_DIRECTORY) + #Space was filling up, so since we've already installed, we can delete the build directory. + echo "rm $(BUILD_DIRECTORY)" + rm -r $(BUILD_DIRECTORY) + echo "ls" ls $(PipExecutable) install numpy pytest echo "$(PythonExecutable) ./setup.py bdist_wheel" $(PythonExecutable) ./setup.py bdist_wheel + echo "rm -r build/" + rm -r build/ + echo "$(PythonExecutable) ./setup_rrplugins.py bdist_wheel" $(PythonExecutable) ./setup_rrplugins.py bdist_wheel # cd dist @@ -1113,12 +1173,11 @@ stages: # echo "wheel is: $wheel" # $(PipExecutable) install $wheel displayName: Build Pip Wheel - - task: CopyFiles@2 - inputs: - sourceFolder: '$(INSTALL_DIRECTORY)' - contents: '**' - targetFolder: '$(Build.ArtifactStagingDirectory)/roadrunner-manylinux2014-$(PythonName)' - displayName: Copy install to artifact staging area + - script: | + echo "mv $(INSTALL_DIRECTORY)/* $(Build.ArtifactStagingDirectory)/roadrunner-manylinux2014-$(PythonName)" + mkdir $(Build.ArtifactStagingDirectory)/roadrunner-manylinux2014-$(PythonName) + mv $(INSTALL_DIRECTORY)/* $(Build.ArtifactStagingDirectory)/roadrunner-manylinux2014-$(PythonName) + displayName: Move install to staging directory - task: PublishBuildArtifacts@1 inputs: pathToPublish: '$(Build.ArtifactStagingDirectory)/roadrunner-manylinux2014-$(PythonName)' diff --git a/docker/INSTRUCTIONS.txt b/docker/INSTRUCTIONS.txt new file mode 100644 index 0000000000..00be8a0e7e --- /dev/null +++ b/docker/INSTRUCTIONS.txt @@ -0,0 +1,22 @@ +When updating (i.e. when Python comes out with a new version): + +* Change roadrunner-manylinux2014/Dockerfile to reflect new python versions. +* Make sure docker is running (launching it on the lab mac works) +* Don't actually use the GUI; just use the command line, once it's running. +* Run the following commands: + +cd roadrunner/docker/ +docker build roadrunner-manylinux2014 +docker images + +* Note the IMAGE ID of the thing you just built (We'll say it's 5555): + +docker tag 5555 sysbiouw/roadrunner-manylinux2014 +docker login -u "sysbiouw" docker.io +docker push sysbiouw/roadrunner-manylinux2014 + + + +If you need to update one of the other images, follow the same instructions, but build up from -base to -add-deps, depending on what needs to be updated. (i.e. if the deps need to be updated, do that first, then the Python with the one after that) + +Make sure you build and push roadrunner-manylinux2014-base first, then roadrunner-manylinux2014-add-deps diff --git a/docker/roadrunner-manylinux2014-base/Dockerfile b/docker/roadrunner-manylinux2014-base/Dockerfile index c1691289b6..3c63acc500 100644 --- a/docker/roadrunner-manylinux2014-base/Dockerfile +++ b/docker/roadrunner-manylinux2014-base/Dockerfile @@ -17,9 +17,9 @@ RUN cd gcc-11.1.0 && ./contrib/download_prerequisites RUN cd gcc-11.1.0 && mkdir build && cd build && ../configure --prefix=/usr/local/gcc-11.1.0 --disable-multilib RUN cd gcc-11.1.0/build && make -j 12 CPPFLAGS=-O3 CXXFLAGS=-O3 RUN cd gcc-11.1.0/build && make install -RUN cd .. && rm -rf gcc-11.1.0 +RUN cd .. && rm -rf gcc-11.1.0 gcc-11.1.0.tar.gz -RUN yum install -y nano +#RUN yum install -y nano RUN yum install -y pcre-devel.x86_64 @@ -29,10 +29,10 @@ RUN mkdir -p $swig_source_dir && wget -q --no-check-certificate https://netcolo RUN tar -xf swig-4.0.2.tar.gz RUN cd $swig_source_dir && ./configure --prefix=$swig_install_dir RUN cd $swig_source_dir && make -j 12 && make install +RUN cd $swig_source_dir && mv $swig_install_dir .. && rm -r -f * && mv ../install-swig* . +RUN rm -r -f swig-4.0.2.tar.gz - -RUN yum install -y ncurses-devel -RUN git config --global core.editor nano && git config --global pager.branch false +#RUN git config --global core.editor nano && git config --global pager.branch false #RUN echo -e "export PATH=\"/usr/local/gcc-11.1.0/bin:$PATH\"" >> ~/.bashrc \ diff --git a/docker/roadrunner-manylinux2014/Dockerfile b/docker/roadrunner-manylinux2014/Dockerfile index 2a71702ce0..5d4877cfbd 100644 --- a/docker/roadrunner-manylinux2014/Dockerfile +++ b/docker/roadrunner-manylinux2014/Dockerfile @@ -11,9 +11,9 @@ FROM sysbiouw/roadrunner-manylinux2014-add-deps RUN wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh RUN bash Miniconda3-latest-Linux-x86_64.sh -b -p /Miniconda3 -RUN /Miniconda3/bin/conda create -y --name py38 python=3.8 pytest RUN /Miniconda3/bin/conda create -y --name py39 python=3.9 pytest RUN /Miniconda3/bin/conda create -y --name py310 python=3.10 pytest RUN /Miniconda3/bin/conda create -y --name py311 python=3.11 pytest -c conda-forge +RUN /Miniconda3/bin/conda create -y --name py312 python=3.12 pytest -c conda-forge RUN /Miniconda3/bin/conda init && bash ~/.bashrc && . ~/.bashrc diff --git a/docs/source/PythonAPIReference/cls_RoadRunner.rst b/docs/source/PythonAPIReference/cls_RoadRunner.rst index 2d2aa112e6..e61ad9559e 100644 --- a/docs/source/PythonAPIReference/cls_RoadRunner.rst +++ b/docs/source/PythonAPIReference/cls_RoadRunner.rst @@ -258,6 +258,17 @@ Selections :rtype: list +.. method:: RoadRunner.[element_id] + :module: RoadRunner + + Access the current value of a model element. All model elements with mathematical meaning (species, compartment, or parameters) can be accessed and set using this method. + >>> rr.S1 + 6.3 + >>> rr.S1 = 2.9 + >>> rr.S1 + 2.9 + + .. method:: RoadRunner.getValue(sel) Returns the value for a given selection. For more information on accepted selection types @@ -268,6 +279,29 @@ Selections :type sel: str or SelectionRecord +.. method:: RoadRunner.setValue(sel, value) + + Sets the value for a given selection. For more information on accepted selection types + see :ref:`Selecting Values `. + + :param str sel: a selection that is either a string or a SelectionRecord that was + obtained from createSelection + :param double value: the value of the selection to be set + + +.. method:: RoadRunner.setValues(keysOrDict, values=None) + + Sets a number of values in the roadrunner object all at once. Use either with the first argument defined as a dictionary, or with both arguments defined, with the first as the keys and the second as the values. + see :ref:`Selecting Values `. + + >>> parameters = {'a': 7, 'b': 8} + >>> rr.setValues(parameters) + >>> rr.setValues(parameters.keys(), parameters.values()) + + :param keysOrDict: either a list of id strings to set, or a dictionary with string keys and numerical values. + :type keysOrDict: list or dict + :param double value: the list of values to use. Must be identical in length to 'keysOrDict', and keysOrDict must not be a dictionary. + .. method:: RoadRunner.getSelectedValues() :module: RoadRunner @@ -306,16 +340,25 @@ Model Access ------------ +.. method:: RoadRunner.getValue(eid) + :module: RoadRunner + + The current value of model elements can be obtained or set by using the ID of the element. This works for all model elements with mathematical meaning (species, compartments, or parameters) + >>> rr.S1 + 6.3 + >>> rr.S1 = 2.9 + >>> rr.S1 + 2.9 + + .. method:: RoadRunner.isModelLoaded() :module: RoadRunner Return True if model was loaded; False otherwise - .. py:attribute:: RoadRunner.model :module: RoadRunner - :annotation: None Get the currently loaded model. The model object contains the entire state of the SBML model. @@ -887,9 +930,9 @@ Easy edit to the model without modifying and reloading sbml files. :param str eid: the ID of the event to be removed - :param bool forceRegenerate: indicate whether the new model is regenerated after this function call + :param bool forceRegenerate: indicate whether the new model is regenerated after this function call. - .. method:: RoadrUNNER.regenerateModel () + .. method:: RoadRunner.regenerateModel () :module RoadRunner Call this method to jit compile any model you've constructed using the modeling editing API. This will make the model ready for simulation. diff --git a/requirements.txt b/requirements.txt index a71d27b534..210ba6aec0 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1 +1 @@ -numpy~=1.24 +numpy~=1.26 diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index 0ef1f61c49..43632ee501 100644 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -225,6 +225,7 @@ if (BUILD_LLVM) llvm/FunctionResolver.cpp llvm/GetEventValuesCodeGen.cpp llvm/GetInitialValuesCodeGen.cpp + llvm/GetPiecewiseTriggerCodeGen.cpp llvm/GetValuesCodeGen.cpp llvm/Jit.cpp llvm/JitFactory.cpp @@ -270,6 +271,7 @@ if (BUILD_LLVM) llvm/GetEventValuesCodeGen.h llvm/GetInitialValueCodeGenBase.h llvm/GetInitialValuesCodeGen.h + llvm/GetPiecewiseTriggerCodeGen.h llvm/GetValueCodeGenBase.h llvm/GetValuesCodeGen.h llvm/Jit.h diff --git a/source/CVODEIntegrator.cpp b/source/CVODEIntegrator.cpp index de36bc284f..b1b29fc3f6 100644 --- a/source/CVODEIntegrator.cpp +++ b/source/CVODEIntegrator.cpp @@ -37,8 +37,7 @@ namespace rr { int cvodeDyDtFcn(realtype t, N_Vector cv_y, N_Vector cv_ydot, void *userData); - int cvodeRootFcn(realtype t, N_Vector y, realtype *gout, void *userData); - + int cvodeEventAndPiecewiseRootFcn(realtype t, N_Vector y, realtype *gout, void *userData); // Sets the value of an element in a N_Vector object inline void SetVector(N_Vector v, int Index, double Value) { @@ -657,11 +656,6 @@ namespace rr { mStateVector = N_VNew_Serial(allocStateVectorSize); variableStepPostEventState.resize(allocStateVectorSize); - -// for (int i = 0; i < allocStateVectorSize; i++) { -// SetVector(mStateVector, i, 0.); -// } -// set mStateVector to the values that are currently in the model auto states = new double[allocStateVectorSize]; mModel->getStateVector(states); @@ -698,12 +692,12 @@ namespace rr { handleCVODEError(err); } - if (mModel->getNumEvents() > 0) { - if ((err = CVodeRootInit(mCVODE_Memory, mModel->getNumEvents(), - cvodeRootFcn)) != CV_SUCCESS) { + if (mModel->getNumEvents() + mModel->getNumPiecewiseTriggers() > 0) { + if ((err = CVodeRootInit(mCVODE_Memory, mModel->getNumEvents() + mModel->getNumPiecewiseTriggers(), + cvodeEventAndPiecewiseRootFcn)) != CV_SUCCESS) { handleCVODEError(err); } - rrLog(Logger::LOG_TRACE) << "CVRootInit executed....."; + rrLog(Logger::LOG_TRACE) << "CVRootInit executed for events....."; } /** @@ -1034,16 +1028,22 @@ namespace rr { // int (*CVRootFn)(realtype t, N_Vector y, realtype *gout, void *user_data) // Cvode calls this to check for event changes - int cvodeRootFcn(realtype time, N_Vector y_vector, realtype *gout, void *user_data) { - CVODEIntegrator *cvInstance = (CVODEIntegrator *) user_data; + int cvodeEventAndPiecewiseRootFcn(realtype time, N_Vector y_vector, realtype* gout, void* user_data) { + CVODEIntegrator* cvInstance = (CVODEIntegrator*)user_data; assert(cvInstance && "user data pointer is NULL on CVODE root callback"); - ExecutableModel *model = cvInstance->mModel; + ExecutableModel* model = cvInstance->mModel; - double *y = NV_DATA_S(y_vector); + double* y = NV_DATA_S(y_vector); - model->getEventRoots(time, y, gout); + if (model->getNumEvents() > 0) { + model->getEventRoots(time, y, gout); + } + + if (model->getNumPiecewiseTriggers() > 0) { + model->getPiecewiseTriggerRoots(time, y, gout + model->getNumEvents()); + } return CV_SUCCESS; } diff --git a/source/CVODEIntegrator.h b/source/CVODEIntegrator.h index 0710520711..ec5dbce489 100644 --- a/source/CVODEIntegrator.h +++ b/source/CVODEIntegrator.h @@ -313,7 +313,7 @@ namespace rr { friend int cvodeDyDtFcn(double t, N_Vector cv_y, N_Vector cv_ydot, void *f_data); - friend int cvodeRootFcn(double t, N_Vector y, double *gout, void *g_data); + friend int cvodeEventAndPiecewiseRootFcn(double t, N_Vector y, double *gout, void *g_data); }; diff --git a/source/GillespieIntegrator.cpp b/source/GillespieIntegrator.cpp index 96cd799e9a..6b8fbb5b35 100644 --- a/source/GillespieIntegrator.cpp +++ b/source/GillespieIntegrator.cpp @@ -11,6 +11,7 @@ #include #include +#include #include #include #include @@ -24,29 +25,33 @@ namespace rr { - static std::uint64_t defaultSeed() - { - Setting seedSetting = Config::getValue(Config::RANDOM_SEED); - std::uint64_t seed; - if (auto int32Val = seedSetting.get_if()) { - seed = (std::uint64_t)*int32Val; - } else if (auto uInt32Val= seedSetting.get_if()){ - seed = (std::uint64_t)*uInt32Val; - } else if (auto int64Val = seedSetting.get_if()){ - seed = (std::uint64_t)*int64Val; - } else if (auto uInt64Val= seedSetting.get_if()){ - seed = *uInt64Val; - } else { - throw std::invalid_argument("GillespieIntegrator::defaultSeed: Seed is of incorrect type."); - } - if (seed < 0) - { - seed = std::chrono::high_resolution_clock::now().time_since_epoch().count(); - } - return seed ; - } - - void GillespieIntegrator::initializeFromModel() { + + static std::uint64_t seedValue(Setting seedSetting) { + std::uint64_t seed; + if (auto int32Val = seedSetting.get_if()) { + seed = (std::uint64_t)*int32Val; + } + else if (auto uInt32Val = seedSetting.get_if()) { + seed = (std::uint64_t)*uInt32Val; + } + else if (auto int64Val = seedSetting.get_if()) { + seed = (std::uint64_t)*int64Val; + } + else if (auto uInt64Val = seedSetting.get_if()) { + seed = *uInt64Val; + } + else { + throw std::invalid_argument("seedSetting is of incorrect type."); + } + + return seed; + } + + static std::uint64_t defaultSeed() { + return seedValue(Config::getValue(Config::RANDOM_SEED)); + } + + void GillespieIntegrator::initializeFromModel() { nReactions = mModel->getNumReactions(); reactionRates = new double[nReactions]; reactionRatesBuffer = new double[nReactions]; @@ -61,29 +66,29 @@ namespace rr assert(floatingSpeciesStart >= 0); - setEngineSeed(getValue("seed").get()); - } + setEngineSeed(getValue("seed")); + } - GillespieIntegrator::GillespieIntegrator(ExecutableModel* m) - : Integrator(m), - timeScale(1.0), - stoichScale(1.0), - reactionRates(nullptr), - reactionRatesBuffer(nullptr), - stateVector(nullptr), - stateVectorRate(nullptr) - { - GillespieIntegrator::resetSettings(); + GillespieIntegrator::GillespieIntegrator(ExecutableModel* m) + : Integrator(m), + timeScale(1.0), + stoichScale(1.0), + reactionRates(nullptr), + reactionRatesBuffer(nullptr), + stateVector(nullptr), + stateVectorRate(nullptr) + { + GillespieIntegrator::resetSettings(); if (mModel) initializeFromModel(); - } + } - GillespieIntegrator::~GillespieIntegrator() - { - if (mModel) { - // like in the constructor, we only delete these components - // if the model was present on construction. + GillespieIntegrator::~GillespieIntegrator() + { + if (mModel) { + // like in the constructor, we only delete these components + // if the model was present on construction. delete[] reactionRates; delete[] reactionRatesBuffer; delete[] stateVector; @@ -93,7 +98,7 @@ namespace rr stateVector = nullptr; stateVectorRate = nullptr; } - } + } void GillespieIntegrator::syncWithModel(ExecutableModel* m) { @@ -135,36 +140,35 @@ namespace rr return "Gillespie Direct Method SSA"; } - Integrator::IntegrationMethod GillespieIntegrator::getIntegrationMethod() const - { - return Integrator::Deterministic; - } - - void GillespieIntegrator::setValue(const std::string& key, Setting val) - { - Integrator::setValue(key, val); - - /* In addition to typically value-setting behavior, some settings require further changes - within CVODE. */ - if (key == "seed") - { - try - { - auto seed = val.getAs(); - setEngineSeed(seed); - } - catch (std::exception& e) - { - std::stringstream ss; - ss << "Could not convert the value \"" << val.get(); - ss << "\" to an unsigned long integer. " << std::endl; - ss << "The seed must be a number between 0 and "; - ss << std::numeric_limits::max(); - ss << "; "<(); + ss << "\" to an unsigned long integer. " << std::endl; + ss << "The seed must be a number between 0 and "; + ss << std::numeric_limits::max(); + ss << "; " << std::endl << "error message: " << e.what() << "."; + throw std::invalid_argument(ss.str()); + } + } + } void GillespieIntegrator::resetSettings() { @@ -178,193 +182,193 @@ namespace rr addSetting("maximum_time_step", Setting(0.0), "Maximum Time Step", "Specifies the maximum absolute value of step size allowed. (double)", "(double) The maximum absolute value of step size allowed."); addSetting("nonnegative", Setting(false), "Non-negative species only", "Prevents species amounts from going negative during a simulation. (bool)", "(bool) Enforce non-negative species constraint."); addSetting("max_output_rows", Setting(Config::getInt(Config::MAX_OUTPUT_ROWS)), "Maximum Output Rows", "For variable step size simulations, the maximum number of output rows produced (int).", "(int) This will set the maximum number of output rows for variable step size integration. This may truncate some simulations that may not reach the desired end time, but prevents massive output for simulations where the variable step size ends up decreasing too much. This setting is ignored when the variable_step_size is false, and is also ignored when the output is being written directly to a file."); - addSetting("maximum_num_steps", Setting(0), "Maximum Number of Steps", - "Specifies the maximum number of steps to be taken by the Gillespie solver before reaching the next reporting time. (int)", - "(int) Maximum number of steps to be taken by the Gillespie solver before reaching the next reporting time."); - } - - double GillespieIntegrator::integrate(double t, double hstep) - { - double tf; - bool singleStep; - bool varStep = getValue("variable_step_size").get(); - auto minTimeStep = getValue("minimum_time_step").get(); - auto maxTimeStep = getValue("maximum_time_step").get(); - int maxNumSteps = getValue("maximum_num_steps").get(); - if (maxTimeStep > minTimeStep) - { - hstep = std::min(hstep, maxTimeStep); - } - - if (varStep) - { - if (minTimeStep > 0.0) - { - if (minTimeStep > hstep) // no step bigger than hstep - { - minTimeStep = hstep; - } - tf = t + minTimeStep; - singleStep = false; - } - else - { - tf = t + hstep; - singleStep = true; - } - } - else - { - tf = t + hstep; - singleStep = false; - } - - assert(hstep > 0 || singleStep && "hstep must be > 0 or we must be taking a single step."); - - rrLog(Logger::LOG_DEBUG) << "ssa(" << t << ", " << tf << ")"; - - // get the initial state std::vector - mModel->setTime(t); - mModel->getStateVector(stateVector); - int step = 0; - - while (singleStep || (t < tf && (maxNumSteps == 0 || step < maxNumSteps))) - { - step++; - // random uniform numbers - double r1 = urand(); - double r2 = urand(); - - assert(r1 > 0 && r1 <= 1 && r2 >= 0 && r2 <= 1); - - // sum of propensities - double s = 0; - - // next time - double tau; - - // get the 'propensity' -- reaction rates - mModel->getReactionRates(nReactions, nullptr, reactionRates); - - // sum the propensity - for (int k = 0; k < nReactions; k++) - { - rrLog(Logger::LOG_DEBUG) << "reac rate: " << k << ": " - << reactionRates[k]; - - // if reaction rate is negative, that means reaction goes in reverse, - // this is fine, we just have to reverse the stoichiometry, - // but still need to sum the absolute value of the propensities - // to get tau. - s += std::abs(reactionRates[k]); - } - - // sample tau - if (s > 0) - { - tau = -log(r1) / s; - } - else - { - // no reaction occurs - return tf; - } - if (!singleStep && t + tau > tf) // if time exhausted, don't allow reaction to proceed - { - return tf; - } - - t = t + tau; - - // select reaction - int reaction = -1; - double sp = 0.0; - - r2 = r2 * s; - for (int i = 0; i < nReactions; ++i) - { - sp += std::abs(reactionRates[i]); - if (r2 < sp) - { - reaction = i; - break; - } - } - - assert(reaction >= 0 && reaction < nReactions); - - // update chemical species - // if rate is negative, means reaction goes in reverse, so - // multiply by sign - double sign = (reactionRates[reaction] > 0) - - (reactionRates[reaction] < 0); - - bool skip = false; - - if ((bool)getValue("nonnegative")) { - // skip reactions which cause species amts to become negative - for (int i = floatingSpeciesStart; i < stateVectorSize; ++i) { - if (stateVector[i] - + getStoich(i - floatingSpeciesStart, reaction) - * stoichScale * sign < 0.0) { - skip = true; - break; - } - } - } - - if (!skip) { - for (int i = floatingSpeciesStart; i < stateVectorSize; ++i) - { - stateVector[i] = stateVector[i] - + getStoich(i - floatingSpeciesStart, reaction) - * stoichScale * sign; - - if (stateVector[i] < 0.0) { - rrLog(Logger::LOG_WARNING) << "Error, negative value of " - << stateVector[i] - << " encountred for floating species " - << mModel->getFloatingSpeciesId(i - floatingSpeciesStart); - t = std::numeric_limits::infinity(); - } - } - } - - // rates could be time dependent - mModel->setTime(t); - mModel->setStateVector(stateVector); - - // events - bool triggered = false; - - mModel->getEventTriggers(eventStatus.size(), nullptr, !eventStatus.empty() ? &eventStatus[0] : nullptr); - for(unsigned char eventStatu : eventStatus) { - if (eventStatu) - triggered = true; - } - - if (triggered) { - applyEvents(t, previousEventStatus); - } - - if (!eventStatus.empty()) - memcpy(&previousEventStatus[0], &eventStatus[0], eventStatus.size()*sizeof(unsigned char)); - - - if (singleStep) - { - return t; - } - } - if (t < tf && maxNumSteps > 0 && step >= maxNumSteps) - { - std::stringstream msg; - msg << "GillespieIntegrator::integrate failed: max number of steps (" - << maxNumSteps << ") reached before desired output at time " << tf << "."; - throw std::runtime_error(msg.str()); - } - return t; - } + addSetting("maximum_num_steps", Setting(0), "Maximum Number of Steps", + "Specifies the maximum number of steps to be taken by the Gillespie solver before reaching the next reporting time. (int)", + "(int) Maximum number of steps to be taken by the Gillespie solver before reaching the next reporting time."); + } + + double GillespieIntegrator::integrate(double t, double hstep) + { + double tf; + bool singleStep; + bool varStep = getValue("variable_step_size").get(); + auto minTimeStep = getValue("minimum_time_step").get(); + auto maxTimeStep = getValue("maximum_time_step").get(); + int maxNumSteps = getValue("maximum_num_steps").get(); + if (maxTimeStep > minTimeStep) + { + hstep = std::min(hstep, maxTimeStep); + } + + if (varStep) + { + if (minTimeStep > 0.0) + { + if (minTimeStep > hstep) // no step bigger than hstep + { + minTimeStep = hstep; + } + tf = t + minTimeStep; + singleStep = false; + } + else + { + tf = t + hstep; + singleStep = true; + } + } + else + { + tf = t + hstep; + singleStep = false; + } + + assert(hstep > 0 || singleStep && "hstep must be > 0 or we must be taking a single step."); + + rrLog(Logger::LOG_DEBUG) << "ssa(" << t << ", " << tf << ")"; + + // get the initial state std::vector + mModel->setTime(t); + mModel->getStateVector(stateVector); + int step = 0; + + while (singleStep || (t < tf && (maxNumSteps == 0 || step < maxNumSteps))) + { + step++; + // random uniform numbers + double r1 = urand(); + double r2 = urand(); + + assert(r1 > 0 && r1 <= 1 && r2 >= 0 && r2 <= 1); + + // sum of propensities + double s = 0; + + // next time + double tau; + + // get the 'propensity' -- reaction rates + mModel->getReactionRates(nReactions, nullptr, reactionRates); + + // sum the propensity + for (int k = 0; k < nReactions; k++) + { + rrLog(Logger::LOG_DEBUG) << "reac rate: " << k << ": " + << reactionRates[k]; + + // if reaction rate is negative, that means reaction goes in reverse, + // this is fine, we just have to reverse the stoichiometry, + // but still need to sum the absolute value of the propensities + // to get tau. + s += std::abs(reactionRates[k]); + } + + // sample tau + if (s > 0) + { + tau = -log(r1) / s; + } + else + { + // no reaction occurs + return tf; + } + if (!singleStep && t + tau > tf) // if time exhausted, don't allow reaction to proceed + { + return tf; + } + + t = t + tau; + + // select reaction + int reaction = -1; + double sp = 0.0; + + r2 = r2 * s; + for (int i = 0; i < nReactions; ++i) + { + sp += std::abs(reactionRates[i]); + if (r2 < sp) + { + reaction = i; + break; + } + } + + assert(reaction >= 0 && reaction < nReactions); + + // update chemical species + // if rate is negative, means reaction goes in reverse, so + // multiply by sign + double sign = (reactionRates[reaction] > 0) + - (reactionRates[reaction] < 0); + + bool skip = false; + + if ((bool)getValue("nonnegative")) { + // skip reactions which cause species amts to become negative + for (int i = floatingSpeciesStart; i < stateVectorSize; ++i) { + if (stateVector[i] + + getStoich(i - floatingSpeciesStart, reaction) + * stoichScale * sign < 0.0) { + skip = true; + break; + } + } + } + + if (!skip) { + for (int i = floatingSpeciesStart; i < stateVectorSize; ++i) + { + stateVector[i] = stateVector[i] + + getStoich(i - floatingSpeciesStart, reaction) + * stoichScale * sign; + + if (stateVector[i] < 0.0) { + rrLog(Logger::LOG_WARNING) << "Error, negative value of " + << stateVector[i] + << " encountred for floating species " + << mModel->getFloatingSpeciesId(i - floatingSpeciesStart); + t = std::numeric_limits::infinity(); + } + } + } + + // rates could be time dependent + mModel->setTime(t); + mModel->setStateVector(stateVector); + + // events + bool triggered = false; + + mModel->getEventTriggers(eventStatus.size(), nullptr, !eventStatus.empty() ? &eventStatus[0] : nullptr); + for (unsigned char eventStatu : eventStatus) { + if (eventStatu) + triggered = true; + } + + if (triggered) { + applyEvents(t, previousEventStatus); + } + + if (!eventStatus.empty()) + memcpy(&previousEventStatus[0], &eventStatus[0], eventStatus.size() * sizeof(unsigned char)); + + + if (singleStep) + { + return t; + } + } + if (t < tf && maxNumSteps > 0 && step >= maxNumSteps) + { + std::stringstream msg; + msg << "GillespieIntegrator::integrate failed: max number of steps (" + << maxNumSteps << ") reached before desired output at time " << tf << "."; + throw std::runtime_error(msg.str()); + } + return t; + } void GillespieIntegrator::testRootsAtInitialTime() { @@ -373,13 +377,13 @@ namespace rr applyEvents(mIntegrationStartTime, initialEventStatus); } - void GillespieIntegrator::applyEvents(double timeEnd, std::vector &prevEventStatus) + void GillespieIntegrator::applyEvents(double timeEnd, std::vector& prevEventStatus) { mModel->applyEvents(timeEnd, prevEventStatus.empty() ? nullptr : &prevEventStatus[0], stateVector, stateVector); } - void GillespieIntegrator::restart(double t0) - { + void GillespieIntegrator::restart(double t0) + { if (!mModel) { return; } @@ -398,31 +402,41 @@ namespace rr { mModel->getStateVector(stateVector); } - } + } - void GillespieIntegrator::setListener(IntegratorListenerPtr) - { - } + void GillespieIntegrator::setListener(IntegratorListenerPtr) + { + } - IntegratorListenerPtr GillespieIntegrator::getListener() - { - return IntegratorListenerPtr(); - } + IntegratorListenerPtr GillespieIntegrator::getListener() + { + return IntegratorListenerPtr(); + } - double GillespieIntegrator::urand() - { - return (double)engine() / (double)std::mt19937::max(); - } + double GillespieIntegrator::urand() + { + return (double)engine() / (double)std::mt19937::max(); + } - void GillespieIntegrator::setEngineSeed(std::uint64_t seed) - { - rrLog(Logger::LOG_INFORMATION) << "Using user specified seed value: " << seed; + void GillespieIntegrator::setEngineSeed(Setting seedSetting) + { + std::uint64_t seed = seedValue(seedSetting); + rrLog(Logger::LOG_INFORMATION) << "Using user specified seed value: " << seed; - // MSVC needs an explicit cast, fail to compile otherwise. - engine.seed((std::int64_t)seed); - } + // Checks if seed is not equal to -1 (the value which is considered as the random seed value) + if (seed != -1) { + engine.seed(seed); + // MSVC needs an explicit cast, fail to compile otherwise. + } + else { + seed = (std::uint64_t)getMicroSeconds(); + engine.seed(seed); + Integrator::setValue("seed", seed); + rrLog(Logger::LOG_INFORMATION) << "Using seed value from the clock: " << seed; + } + } - Solver *GillespieIntegrator::construct(ExecutableModel *executableModel) const { + Solver* GillespieIntegrator::construct(ExecutableModel* executableModel) const { return new GillespieIntegrator(executableModel); } diff --git a/source/GillespieIntegrator.h b/source/GillespieIntegrator.h index 987916a22f..8296b74a1f 100644 --- a/source/GillespieIntegrator.h +++ b/source/GillespieIntegrator.h @@ -140,7 +140,7 @@ namespace rr void applyEvents(double timeEnd, std::vector &prevEventStatus); double urand(); - void setEngineSeed(std::uint64_t seed); + void setEngineSeed(Setting seedSetting); inline double getStoich(uint species, uint reaction) { diff --git a/source/Setting.h b/source/Setting.h index 13ae1cd809..3e77b78aef 100644 --- a/source/Setting.h +++ b/source/Setting.h @@ -214,42 +214,88 @@ namespace rr { const std::type_info &inf = typeInfo(); return visit([&](auto &&val) { if constexpr (std::is_convertible_v) { - // if we have an int - if (auto intValue = std::get_if(&value_)) { - // and its negative, - if (*intValue < 0) { - // we cannot convert to unsigned int or unsigned long - if (typeid(As) == typeid(unsigned int) || typeid(As) == typeid(unsigned long)) { - throw std::bad_variant_access{}; - return As{}; - } + std::ostringstream os; + os << "Cannot retrieve setting value: you have requested the value as a "; + os << "\"" << typeid(As).name() << "\", but the value of the setting is "; + std::ostringstream oss_val; + + // if we have some form of negative signed integer + bool isNegInt = false; + uint64_t uint64val = 0; + int64_t int64val = 0; + if (auto intValue = std::get_if(&value_)) + { + if (*intValue < 0) + { + isNegInt = true; + int64val = *intValue; } + uint64val = *intValue; + oss_val << "\"" << *intValue << "\", which is "; } - - // if we have a long, - if (auto lValue = std::get_if(&value_)) { - // and its negative, we cannot convert to unsigned int or unsigned long - if (*lValue < 0) { - if (typeid(As) == typeid(unsigned int) || typeid(As) == typeid(unsigned long)) { - throw std::bad_variant_access{}; - return As{}; // must be present - } + else if (auto intValue = std::get_if(&value_)) + { + if (*intValue < 0) + { + isNegInt = true; + int64val = *intValue; } - // furthermore, if we have a long which has a value greater than - // that of int32 maximum, we have a problem and throw - if (*lValue > ((std::int64_t) (std::numeric_limits::max)())) { - throw std::bad_variant_access{}; // has annoying private constructor so we can't add arguments - return As{}; // must be present + oss_val << "\"" << *intValue << "\", which is "; + uint64val = *intValue; + } + else if (auto intValue = std::get_if(&value_)) + { + oss_val << "\"" << *intValue << "\", which is "; + uint64val = *intValue; + } + else if (auto intValue = std::get_if(&value_)) + { + oss_val << "\"" << *intValue << "\", which is "; + uint64val = *intValue; + } + if (isNegInt) + { + // we cannot convert to an unsigned type + if (typeid(As) == typeid(unsigned int) || typeid(As) == typeid(unsigned long)) { + // would prefer to throw std::bad_variant_access, but it seems + // it does not have the appropriate constructor (?) + os << oss_val.str() << "negative." << std::endl; + throw std::invalid_argument(os.str()); + return As{}; // oddly enough, this *is* necessary } } // if we have a double, with value greater than std::numeric_limits::max // and we try to convert to float, we have an error if (auto lValue = std::get_if(&value_)) { if (*lValue > ((std::numeric_limits::max)())) { - throw std::bad_variant_access{}; + os << "\"" << *lValue << "\", which is too large." << std::endl; + throw std::invalid_argument(os.str()); return As{}; // must be present } } + //Now check if the value is larger than the max value for the type: + if (typeid(As) == typeid(int) && + (uint64val > ((std::int64_t)(std::numeric_limits::max)()) && !isNegInt) || + (isNegInt && int64val < ((std::int64_t)(std::numeric_limits::min)()))) + { + os << oss_val.str() << "too large." << std::endl; + throw std::invalid_argument(os.str()); + return As{}; // must be present + } + else if (typeid(As) == typeid(unsigned int) && + uint64val > ((std::int64_t)(std::numeric_limits::max)())) + { + os << oss_val.str() << "too large." << std::endl; + throw std::invalid_argument(os.str()); + return As{}; // must be present + } + else if (typeid(As) == typeid(int64_t) && + (uint64val > ((std::int64_t)(std::numeric_limits::max)()) && !isNegInt)) + { + os << oss_val.str() << "too large." << std::endl; + throw std::invalid_argument(os.str()); + return As{}; // must be present + } return As(val); } else { std::ostringstream os; diff --git a/source/llvm/ASTNodeCodeGen.cpp b/source/llvm/ASTNodeCodeGen.cpp index ca7ea366f1..b2f4174c31 100644 --- a/source/llvm/ASTNodeCodeGen.cpp +++ b/source/llvm/ASTNodeCodeGen.cpp @@ -241,7 +241,9 @@ llvm::Value* ASTNodeCodeGen::codeGen(const libsbml::ASTNode* ast) case AST_FUNCTION_DELAY: result = delayExprCodeGen(ast); break; - + case AST_FUNCTION_RATE_OF: + result = rateOfCodeGen(ast); + break; case AST_LAMBDA: result = notImplemented(ast); break; @@ -637,6 +639,159 @@ llvm::Value* ASTNodeCodeGen::delayExprCodeGen(const libsbml::ASTNode* ast) //return codeGen(ast->getChild(0)); } +llvm::Value* ASTNodeCodeGen::rateOfCodeGen(const libsbml::ASTNode* ast) +{ + if (ast->getNumChildren() != 1) { + throw_llvm_exception("AST type 'rateOf' requires exactly one child."); + } + + ASTNode* child = ast->getChild(0); + if (child->getType() != AST_NAME) + { + char* formula = SBML_formulaToL3String(ast); + std::stringstream err; + err << "The rateOf csymbol may only be used on individual symbols, i.e. 'rateOf(S1)'. The expression '" << formula << "' is illegal."; + free(formula); + throw_llvm_exception(err.str()) + } + string name = child->getName(); + const LLVMModelDataSymbols& dataSymbols = ctx.getModelDataSymbols(); + //ModelDataIRBuilder mdbuilder(modelData, dataSymbols, builder); + Value* rate = NULL; + resolver.recursiveSymbolPush("rateOf(" + name + ")"); + //Looking for a rate. Our options are: local parameter, floating species, rate rule target, assignment rule target (an error), and everything else is zero. + if (resolver.isLocalParameter(name)) + { + rate = ConstantFP::get(builder.getContext(), APFloat(0.0)); + } + else if (dataSymbols.isIndependentFloatingSpecies(name)) + { + //NOTE: if we stored rates persistently, we could use the following instead: + //rate = mdbuilder.createFloatSpeciesAmtRateLoad(name, name + "_amtRate"); + const Model* model = ctx.getModel(); + const Species* species = model->getSpecies(name); + bool speciesIsAmount = species->getHasOnlySubstanceUnits(); + string compId = species->getCompartment(); + bool compIsConst = true; + + if (!speciesIsAmount) { + //Determine if the compartment changes in time. + const Rule* compRule = model->getRule(compId); + if (compRule != NULL) { + if (compRule->getTypeCode() == SBML_ASSIGNMENT_RULE) { + std::stringstream err; + err << "Unable to calculate rateOf(" << name << "), because " + << name << " is a concentration, not an amount, but its compartment (" + << compId << ") chages due to an assignment rule. Since we cannot" + << " calculate derivatives on the fly, we cannot determine the rate of" + << " change of the species " << name << " concentration"; + throw_llvm_exception(err.str()); + } + else { + //The only other option is a rate rule, which means the compartment changes in time. + compIsConst = false; + } + } + } + + //The functions are different depending on how complicated things get, so we need references. The final Value will be created from the overallRef ASTNode. + ASTNode amtRate(AST_PLUS); + ASTNode* amtRateRef = &amtRate; + ASTNode* overallRef = &amtRate; + + //Deal with conversion factors: + string conversion_factor = ""; + if (model->isSetConversionFactor()) { + conversion_factor = model->getConversionFactor(); + } + if (species->isSetConversionFactor()) { + conversion_factor = species->getConversionFactor(); + } + if (!conversion_factor.empty()) { + amtRate.setType(AST_TIMES); //rate * conversion factor + ASTNode* unscaledRate = new ASTNode(AST_PLUS); + amtRate.addChild(unscaledRate); + ASTNode* cf = new ASTNode(AST_NAME); + cf->setName(conversion_factor.c_str()); + amtRate.addChild(cf); + amtRateRef = unscaledRate; + } + + //Sum all the kinetic laws * stoichiometries in which this species appears: + for (int rxn = 0; rxn < model->getNumReactions(); rxn++) { + const Reaction* reaction = model->getReaction(rxn); + if (reaction->getReactant(name) == NULL && + reaction->getProduct(name) == NULL) { + continue; + } + ASTNode* times = new ASTNode(AST_TIMES); + ASTNode* stoich = new ASTNode(AST_NAME); + stoich->setName((reaction->getId() + ":" + name).c_str()); + times->addChild(stoich); + times->addChild(reaction->getKineticLaw()->getMath()->deepCopy()); + amtRateRef->addChild(times); + } + + //Now deal with the case where the species symbol is a concentration, not an amount. + ASTNode scaledConcentrationRate(AST_DIVIDE); + ASTNode combinedConcentrationRate(AST_MINUS); + if (!speciesIsAmount) { + //Need to divide by the compartment size (rate / comp) + scaledConcentrationRate.addChild(overallRef->deepCopy()); + ASTNode* compsize = new ASTNode(AST_NAME); + compsize->setName(compId.c_str()); + scaledConcentrationRate.addChild(compsize); + overallRef = &scaledConcentrationRate; + + //If the compartment changes in time, it becomes even more complicated! + if (!compIsConst) { + //The equation is: rateOf([S1]) = rateOf(S1)/C - [S1]/C * rateOf(C) + combinedConcentrationRate.addChild(overallRef->deepCopy()); + ASTNode* subdiv = new ASTNode(AST_DIVIDE); + ASTNode* mult = new ASTNode(AST_TIMES); + ASTNode* species = new ASTNode(AST_NAME); + species->setName(name.c_str()); + ASTNode* rateOfComp = new ASTNode(AST_FUNCTION_RATE_OF); + ASTNode* comp = new ASTNode(AST_NAME); + comp->setName(compId.c_str()); + rateOfComp->addChild(comp); + mult->addChild(rateOfComp); + mult->addChild(species); + subdiv->addChild(mult); + subdiv->addChild(comp->deepCopy()); + combinedConcentrationRate.addChild(subdiv); + overallRef = &combinedConcentrationRate; + } + + } + //string formula = SBML_formulaToL3String(overallRef); + rate = ASTNodeCodeGen(builder, resolver, ctx, modelData).codeGenDouble(overallRef); + } + else if (dataSymbols.hasRateRule(name)) + { + //NOTE: if we stored rates persistently, we could use the following instead: + //rate = mdbuilder.createRateRuleRateLoad(name, name + "_rate"); + const LLVMModelSymbols& modelSymbols = ctx.getModelSymbols(); + SymbolForest::ConstIterator i = modelSymbols.getRateRules().find( + name); + if (i != modelSymbols.getRateRules().end()) + { + rate = ASTNodeCodeGen(builder, resolver, ctx, modelData).codeGenDouble(i->second); + } + } + else if (dataSymbols.hasAssignmentRule(name)) + { + throw_llvm_exception("Unable to define the rateOf for symbol '" + name + "' as it is changed by an assignment rule."); + } + else + { + rate = ConstantFP::get(builder.getContext(), APFloat(0.0)); + } + assert(rate); + resolver.recursiveSymbolPop(); + return rate; +} + llvm::Value* ASTNodeCodeGen::nameExprCodeGen(const libsbml::ASTNode* ast) { switch(ast->getType()) diff --git a/source/llvm/ASTNodeCodeGen.h b/source/llvm/ASTNodeCodeGen.h index c78cd9cdd0..5153ca0ff7 100644 --- a/source/llvm/ASTNodeCodeGen.h +++ b/source/llvm/ASTNodeCodeGen.h @@ -51,6 +51,8 @@ class ASTNodeCodeGen llvm::Value *delayExprCodeGen(const libsbml::ASTNode *ast); + llvm::Value* rateOfCodeGen(const libsbml::ASTNode* ast); + llvm::Value *nameExprCodeGen(const libsbml::ASTNode *ast); llvm::Value *realExprCodeGen(const libsbml::ASTNode *ast); diff --git a/source/llvm/CodeGen.h b/source/llvm/CodeGen.h index 4bdac66495..caf5c55083 100644 --- a/source/llvm/CodeGen.h +++ b/source/llvm/CodeGen.h @@ -50,6 +50,8 @@ namespace rrllvm virtual void recursiveSymbolPop() = 0; + virtual bool isLocalParameter(const std::string& symbol) = 0; + /** * nested conditionals (or functions?) can push a local cache block, where * symbols would be chached. These need to be popped as these symbols are diff --git a/source/llvm/CodeGenBase.h b/source/llvm/CodeGenBase.h index e11c62f0a9..dde40640c6 100644 --- a/source/llvm/CodeGenBase.h +++ b/source/llvm/CodeGenBase.h @@ -78,7 +78,6 @@ class CodeGenBase builder(*mgc.getJitNonOwning()->getBuilderNonOwning()), options(mgc.getOptions()), function(0) -// functionPassManager(mgc.getJitNonOwning().getFunctionPassManager()) { }; diff --git a/source/llvm/EvalVolatileStoichCodeGen.cpp b/source/llvm/EvalVolatileStoichCodeGen.cpp index d2f3989419..4dadfb94ab 100644 --- a/source/llvm/EvalVolatileStoichCodeGen.cpp +++ b/source/llvm/EvalVolatileStoichCodeGen.cpp @@ -42,6 +42,8 @@ namespace rrllvm { ModelDataLoadSymbolResolver resolver(modelData, modelGenContext); ModelDataIRBuilder mdbuilder(modelData, dataSymbols, builder); + LLVMModelDataSymbols* dataSymbolsPtr = const_cast(&dataSymbols); + ASTNodeCodeGen astCodeGen(builder, resolver, modelGenContext, modelData); const ListOfReactions *reactions = model->getListOfReactions(); @@ -79,7 +81,7 @@ namespace rrllvm { assert(value && "value for species reference stoichiometry is 0"); const LLVMModelDataSymbols::SpeciesReferenceInfo &info = - dataSymbols.getNamedSpeciesReferenceInfo(p->getId()); + dataSymbolsPtr->getNamedSpeciesReferenceInfo(p->getId()); mdbuilder.createStoichiometryStore(info.row, info.column, value, p->getId()); @@ -117,7 +119,7 @@ namespace rrllvm { value = builder.CreateFMul(negOne, value, "neg_" + r->getId()); const LLVMModelDataSymbols::SpeciesReferenceInfo &info = - dataSymbols.getNamedSpeciesReferenceInfo(r->getId()); + dataSymbolsPtr->getNamedSpeciesReferenceInfo(r->getId()); mdbuilder.createStoichiometryStore(info.row, info.column, value, r->getId()); diff --git a/source/llvm/FunctionResolver.cpp b/source/llvm/FunctionResolver.cpp index e24a918515..378bf96c6c 100644 --- a/source/llvm/FunctionResolver.cpp +++ b/source/llvm/FunctionResolver.cpp @@ -126,4 +126,9 @@ void FunctionResolver::recursiveSymbolPop() parentResolver.recursiveSymbolPop(); } +bool FunctionResolver::isLocalParameter(const std::string& symbol) +{ + return false; +} + } /* namespace rr */ diff --git a/source/llvm/FunctionResolver.h b/source/llvm/FunctionResolver.h index e22d59fcf6..44a53f47a4 100644 --- a/source/llvm/FunctionResolver.h +++ b/source/llvm/FunctionResolver.h @@ -35,6 +35,8 @@ class FunctionResolver: public LoadSymbolResolver virtual void recursiveSymbolPop(); + virtual bool isLocalParameter(const std::string& symbol); + private: LoadSymbolResolver& parentResolver; const ModelGeneratorContext& modelGenContext; diff --git a/source/llvm/GetPiecewiseTriggerCodeGen.cpp b/source/llvm/GetPiecewiseTriggerCodeGen.cpp new file mode 100644 index 0000000000..b8f8b6f390 --- /dev/null +++ b/source/llvm/GetPiecewiseTriggerCodeGen.cpp @@ -0,0 +1,110 @@ +/* + * GetPiecewiseTriggerCodeGen.cpp + * + * Created on: Aug 10, 2013 + * Author: andy + */ +#pragma hdrstop +#include "GetPiecewiseTriggerCodeGen.h" +#include "LLVMException.h" +#include "ModelDataSymbolResolver.h" +#include "rrLogger.h" + +#include +#include +#include + +using namespace llvm; + +using namespace libsbml; + +namespace rrllvm +{ + const char* GetPiecewiseTriggerCodeGen::FunctionName = "getPiecewiseTrigger"; + const char* GetPiecewiseTriggerCodeGen::IndexArgName = "triggerIndx"; + + GetPiecewiseTriggerCodeGen::GetPiecewiseTriggerCodeGen(const ModelGeneratorContext& mgc) + : CodeGenBase(mgc) + , piecewiseTriggers(mgc.getPiecewiseTriggers()) + { + }; + + llvm::Value* GetPiecewiseTriggerCodeGen::codeGen() + { + // make the set init value function + llvm::Type* argTypes[] = { + llvm::PointerType::get(ModelDataIRBuilder::getStructType(this->module), 0), + llvm::Type::getInt32Ty(this->context) + }; + + const char* argNames[] = { + "modelData", IndexArgName + }; + + llvm::Value* args[] = { 0, 0 }; + + llvm::Type* retType = getRetType(); + + llvm::BasicBlock* entry = this->codeGenHeader(FunctionName, retType, + argTypes, argNames, args); + + ModelDataLoadSymbolResolver resolver(args[0], this->modelGenContext); + + ASTNodeCodeGen astCodeGen(this->builder, resolver, this->modelGenContext, args[0]); + + // default, return NaN + llvm::BasicBlock* def = llvm::BasicBlock::Create(this->context, "default", this->function); + this->builder.SetInsertPoint(def); + + llvm::Value* defRet = createRet(0); + this->builder.CreateRet(defRet); + + // write the switch at the func entry point, the switch is also the + // entry block terminator + this->builder.SetInsertPoint(entry); + + llvm::SwitchInst* s = this->builder.CreateSwitch(args[1], def, piecewiseTriggers->size()); + + for (uint i = 0; i < piecewiseTriggers->size(); ++i) + { + char block_name[64]; + std::sprintf(block_name, "piecewiseTrigger_%i_block", i); + llvm::BasicBlock* block = llvm::BasicBlock::Create(this->context, block_name, this->function); + this->builder.SetInsertPoint(block); + resolver.flushCache(); + + llvm::Value* value = astCodeGen.codeGenBoolean((*piecewiseTriggers)[i]); + + // convert type to return type + value = createRet(value); + + this->builder.CreateRet(value); + s->addCase(llvm::ConstantInt::get(llvm::Type::getInt32Ty(this->context), i), block); + } + + return this->verifyFunction(); + } + + llvm::Type* GetPiecewiseTriggerCodeGen::getRetType() + { + return llvm::Type::getInt8Ty(context); + }; + + llvm::Value* GetPiecewiseTriggerCodeGen::createRet(llvm::Value* value) + { + if (!value) + { + return llvm::ConstantInt::get(getRetType(), 0xff, false); + } + else + { + return builder.CreateIntCast(value, getRetType(), false); + } + } + + +} /* namespace rr */ + + + + diff --git a/source/llvm/GetPiecewiseTriggerCodeGen.h b/source/llvm/GetPiecewiseTriggerCodeGen.h new file mode 100644 index 0000000000..4e4cd6e846 --- /dev/null +++ b/source/llvm/GetPiecewiseTriggerCodeGen.h @@ -0,0 +1,59 @@ +/* + * GetPiecewiseTriggerCodeGen.h + * + * Created on: Aug 10, 2013 + * Author: andy + */ + +#ifndef RRLLVMGetPiecewiseTriggerCodeGen_H_ +#define RRLLVMGetPiecewiseTriggerCodeGen_H_ + +#include "CodeGenBase.h" +#include "ModelGeneratorContext.h" +#include "SymbolForest.h" +#include "ASTNodeCodeGen.h" +#include "ASTNodeFactory.h" +#include "ModelDataIRBuilder.h" +#include "ModelDataSymbolResolver.h" +#include "LLVMException.h" +#include "rrLogger.h" +#include +#include +#include +#include + +namespace rrllvm +{ + //Based on GetEventTriggerCodeGen (-LS) + + typedef unsigned char (*GetPiecewiseTriggerCodeGen_FunctionPtr)(LLVMModelData*, size_t); + + class GetPiecewiseTriggerCodeGen : + public CodeGenBase + { + public: + GetPiecewiseTriggerCodeGen(const ModelGeneratorContext& mgc); + virtual ~GetPiecewiseTriggerCodeGen() {}; + + llvm::Value* codeGen(); + + typedef GetPiecewiseTriggerCodeGen_FunctionPtr FunctionPtr; + + static const char* FunctionName; + static const char* IndexArgName; + + llvm::Type* getRetType(); + + llvm::Value* createRet(llvm::Value*); + + private: + const std::vector* piecewiseTriggers; + }; + + +} /* namespace rr */ + + + + +#endif /* RRLLVMGETVALUECODEGENBASE_H_ */ diff --git a/source/llvm/Jit.cpp b/source/llvm/Jit.cpp index d7fd3a883f..72aab7b5e6 100644 --- a/source/llvm/Jit.cpp +++ b/source/llvm/Jit.cpp @@ -128,6 +128,8 @@ namespace rrllvm { "eventTrigger")); modelResources->eventAssignPtr = reinterpret_cast(lookupFunctionAddress( "eventAssign")); + modelResources->getPiecewiseTriggerPtr = reinterpret_cast(lookupFunctionAddress( + "getPiecewiseTrigger")); modelResources->evalVolatileStoichPtr = reinterpret_cast(lookupFunctionAddress( "evalVolatileStoich")); modelResources->evalConversionFactorPtr = reinterpret_cast(lookupFunctionAddress( diff --git a/source/llvm/KineticLawParameterResolver.cpp b/source/llvm/KineticLawParameterResolver.cpp index 7fcf645016..37b20defa8 100644 --- a/source/llvm/KineticLawParameterResolver.cpp +++ b/source/llvm/KineticLawParameterResolver.cpp @@ -67,4 +67,12 @@ void KineticLawParameterResolver::recursiveSymbolPop() parentResolver.recursiveSymbolPop(); } +bool KineticLawParameterResolver::isLocalParameter(const std::string& symbol) +{ + if (kineticLaw.getLocalParameter(symbol) != NULL) { + return true; + } + return kineticLaw.getParameter(symbol) != NULL; +} + } /* namespace rr */ diff --git a/source/llvm/KineticLawParameterResolver.h b/source/llvm/KineticLawParameterResolver.h index 66783d954b..6de83e04aa 100644 --- a/source/llvm/KineticLawParameterResolver.h +++ b/source/llvm/KineticLawParameterResolver.h @@ -33,6 +33,8 @@ namespace rrllvm virtual void recursiveSymbolPop(); + virtual bool isLocalParameter(const std::string& symbol); + private: LoadSymbolResolver& parentResolver; const libsbml::KineticLaw& kineticLaw; diff --git a/source/llvm/LLVMExecutableModel.cpp b/source/llvm/LLVMExecutableModel.cpp index f7b40932b3..229d97d9a5 100644 --- a/source/llvm/LLVMExecutableModel.cpp +++ b/source/llvm/LLVMExecutableModel.cpp @@ -182,6 +182,7 @@ LLVMExecutableModel::LLVMExecutableModel() : getEventDelayPtr(0), eventTriggerPtr(0), eventAssignPtr(0), + getPiecewiseTriggerPtr(0), evalVolatileStoichPtr(0), evalConversionFactorPtr(0), setBoundarySpeciesAmountPtr(0), @@ -229,6 +230,7 @@ LLVMExecutableModel::LLVMExecutableModel( getEventDelayPtr(modelResources->getEventDelayPtr), eventTriggerPtr(modelResources->eventTriggerPtr), eventAssignPtr(modelResources->eventAssignPtr), + getPiecewiseTriggerPtr(modelResources->getPiecewiseTriggerPtr), evalVolatileStoichPtr(modelResources->evalVolatileStoichPtr), evalConversionFactorPtr(modelResources->evalConversionFactorPtr), setBoundarySpeciesAmountPtr(modelResources->setBoundarySpeciesAmountPtr), @@ -294,6 +296,7 @@ LLVMExecutableModel::LLVMExecutableModel(std::istream& in, uint modelGeneratorOp getEventDelayPtr = resources->getEventDelayPtr; eventTriggerPtr = resources->eventTriggerPtr; eventAssignPtr = resources->eventAssignPtr; + getPiecewiseTriggerPtr = resources->getPiecewiseTriggerPtr; evalVolatileStoichPtr = resources->evalVolatileStoichPtr; evalConversionFactorPtr = resources->evalConversionFactorPtr; setBoundarySpeciesAmountPtr = resources->setBoundarySpeciesAmountPtr; @@ -1871,6 +1874,11 @@ void LLVMExecutableModel::getEventIds(std::list& out) std::copy(eventIds.begin(), eventIds.end(), std::back_inserter(out)); } +int LLVMExecutableModel::getNumPiecewiseTriggers() +{ + return modelData->numPiecewiseTriggers; +} + void LLVMExecutableModel::getAssignmentRuleIds(std::list& out) { std::vector arIds = symbols->getAssignmentRuleIds(); @@ -2462,12 +2470,6 @@ void LLVMExecutableModel::getEventRoots(double time, const double* y, double* g if (y) { - //memcpy(modelData->rateRuleValues, y, - // modelData->numRateRules * sizeof(double)); - - //memcpy(modelData->floatingSpeciesAmounts, y + modelData->numRateRules, - // modelData->numIndFloatingSpecies * sizeof(double)); - modelData->rateRuleValuesAlias = const_cast(y); modelData->floatingSpeciesAmountsAlias = const_cast(y + modelData->numRateRules); @@ -2487,6 +2489,34 @@ void LLVMExecutableModel::getEventRoots(double time, const double* y, double* g return; } +void LLVMExecutableModel::getPiecewiseTriggerRoots(double time, const double* y, double* gdot) +{ + modelData->time = time; + + double* savedRateRules = modelData->rateRuleValuesAlias; + double* savedFloatingSpeciesAmounts = modelData->floatingSpeciesAmountsAlias; + + if (y) + { + modelData->rateRuleValuesAlias = const_cast(y); + modelData->floatingSpeciesAmountsAlias = const_cast(y + modelData->numRateRules); + + evalVolatileStoichPtr(modelData); + } + + for (uint i = 0; i < modelData->numPiecewiseTriggers; ++i) + { + unsigned char triggered = getPiecewiseTriggerPtr(modelData, i); + + gdot[i] = triggered ? 1.0 : -1.0; + } + + modelData->rateRuleValuesAlias = savedRateRules; + modelData->floatingSpeciesAmountsAlias = savedFloatingSpeciesAmounts; + + return; +} + double LLVMExecutableModel::getNextPendingEventTime(bool pop) { return pendingEvents.getNextPendingEventTime(); diff --git a/source/llvm/LLVMExecutableModel.h b/source/llvm/LLVMExecutableModel.h index ea82993979..30b4b26cc2 100644 --- a/source/llvm/LLVMExecutableModel.h +++ b/source/llvm/LLVMExecutableModel.h @@ -32,6 +32,7 @@ #include "GetEventValuesCodeGen.h" #include "EventAssignCodeGen.h" #include "EventTriggerCodeGen.h" +#include "GetPiecewiseTriggerCodeGen.h" #include "EvalVolatileStoichCodeGen.h" #include "EvalConversionFactorCodeGen.h" #include "SetValuesCodeGen.h" @@ -97,7 +98,7 @@ class RR_DECLSPEC LLVMExecutableModel: public rr::ExecutableModel * evaluate the initial conditions specified in the sbml, this entails * evaluating all InitialAssigments, AssigmentRules, initial values, etc... * - * The the model state is fully set. + * Then the model state is fully set. */ void evalInitialConditions(uint32_t flags = 0); @@ -489,6 +490,8 @@ class RR_DECLSPEC LLVMExecutableModel: public rr::ExecutableModel virtual void getEventRoots(double time, const double* y, double* gdot); + virtual void getPiecewiseTriggerRoots(double time, const double* y, double* gdot); + virtual double getNextPendingEventTime(bool pop); virtual int getPendingEventSize(); @@ -563,6 +566,9 @@ class RR_DECLSPEC LLVMExecutableModel: public rr::ExecutableModel virtual int getEventIndex(const std::string& eid); virtual std::string getEventId(size_t index); virtual void getEventIds(std::list& out); + + virtual int getNumPiecewiseTriggers(); + virtual void getAssignmentRuleIds(std::list& out); virtual void getRateRuleIds(std::list& out); virtual void getInitialAssignmentIds(std::list& out); @@ -656,6 +662,7 @@ class RR_DECLSPEC LLVMExecutableModel: public rr::ExecutableModel GetEventDelayCodeGen::FunctionPtr getEventDelayPtr; EventTriggerCodeGen::FunctionPtr eventTriggerPtr; EventAssignCodeGen::FunctionPtr eventAssignPtr; + GetPiecewiseTriggerCodeGen::FunctionPtr getPiecewiseTriggerPtr; EvalVolatileStoichCodeGen::FunctionPtr evalVolatileStoichPtr; EvalConversionFactorCodeGen::FunctionPtr evalConversionFactorPtr; diff --git a/source/llvm/LLVMModelData.cpp b/source/llvm/LLVMModelData.cpp index b6fd02e786..73ab8f1a20 100644 --- a/source/llvm/LLVMModelData.cpp +++ b/source/llvm/LLVMModelData.cpp @@ -121,6 +121,7 @@ void LLVMModelData_save(LLVMModelData *data, std::ostream& out) rr::saveBinary(out, data->numInitGlobalParameters); rr::saveBinary(out, data->numEvents); + rr::saveBinary(out, data->numPiecewiseTriggers); rr::saveBinary(out, data->stateVectorSize); //Save the stoichiometry matrix rr::csr_matrix_dump_binary(data->stoichiometry, out); @@ -199,8 +200,9 @@ LLVMModelData* LLVMModelData_from_save(std::istream& in) rr::loadBinary(in, data->numInitGlobalParameters); rr::loadBinary(in, data->numEvents); + rr::loadBinary(in, data->numPiecewiseTriggers); rr::loadBinary(in, data->stateVectorSize); - + //Load the stoichiometry matrix data->stoichiometry = rr::csr_matrix_new_from_binary(in); diff --git a/source/llvm/LLVMModelData.h b/source/llvm/LLVMModelData.h index a7a2a1fa67..0a789a262f 100644 --- a/source/llvm/LLVMModelData.h +++ b/source/llvm/LLVMModelData.h @@ -117,23 +117,26 @@ struct LLVMModelData //Event stuff unsigned numEvents; // 15 + //Piecewise triggers that get treated like events for rootfinding purposes. + unsigned numPiecewiseTriggers; // 16 + /** * number of items in the state std::vector. * should be numIndFloatingSpecies + numRateRules */ - unsigned stateVectorSize; // 16 + unsigned stateVectorSize; // 17 /** * the state std::vector, this is usually a pointer to a block of data * owned by the integrator. */ - double* stateVector; // 17 + double* stateVector; // 18 /** * the rate of change of the state std::vector, this is usually a pointer to * a block of data owned by the integrator. */ - double* stateVectorRate; // 18 + double* stateVectorRate; // 19 /** * the rate of change of all elements who's dynamics are determined @@ -144,7 +147,7 @@ struct LLVMModelData * * Normally NULL, only valid durring an evalModel call. */ - double* rateRuleRates; // 19 + double* rateRuleRates; // 20 @@ -154,7 +157,7 @@ struct LLVMModelData * This pointer is ONLY valid during an evalModel call, otherwise it is * zero. TODO, this needs be be moved to a parameter. */ - double* floatingSpeciesAmountRates; // 20 + double* floatingSpeciesAmountRates; // 21 @@ -170,8 +173,8 @@ struct LLVMModelData * units: volume */ - double* compartmentVolumesAlias; // 21 - double* initCompartmentVolumesAlias; // 22 + double* compartmentVolumesAlias; // 22 + double* initCompartmentVolumesAlias; // 23 /** @@ -182,16 +185,16 @@ struct LLVMModelData * Note that dependent floating species which have a rate rule will not be stored * in this block, instead, they will be stored in RateRule block */ - double* initFloatingSpeciesAmountsAlias; // 23 + double* initFloatingSpeciesAmountsAlias; // 24 - double* boundarySpeciesAmountsAlias; // 24 - double* initBoundarySpeciesAmountsAlias; // 25 + double* boundarySpeciesAmountsAlias; // 25 + double* initBoundarySpeciesAmountsAlias; // 26 - double* globalParametersAlias; // 26 - double* initGlobalParametersAlias; // 27 + double* globalParametersAlias; // 27 + double* initGlobalParametersAlias; // 28 - double* reactionRatesAlias; // 28 + double* reactionRatesAlias; // 29 /** * All of the elelments which have a rate rule are stored here, including @@ -207,7 +210,7 @@ struct LLVMModelData * of this struct. * */ - double* rateRuleValuesAlias; // 29 + double* rateRuleValuesAlias; // 30 @@ -218,22 +221,22 @@ struct LLVMModelData * This pointer is part of the state std::vector. When any function is called by * CVODE, this is actually a pointer to a CVODE owned memory block. */ - double* floatingSpeciesAmountsAlias; // 30 + double* floatingSpeciesAmountsAlias; // 31 /** * binary data layout: * - * compartmentVolumes [numIndCompartmentVolumes] // 31 - * initCompartmentVolumes [numInitCompartmentVolumes] // 32 - * initFloatingSpeciesAmounts [numInitFloatingSpecies] // 33 - * boundarySpeciesAmounts [numIndBoundarySpecies] // 34 - * initBoundarySpeciesAmounts [numInitBoundarySpecies] // 35 - * globalParameters [numIndGlobalParameters] // 36 - * initGlobalParameters [numInitGlobalParameters] // 37 - * reactionRates [numReactions] // 38 + * compartmentVolumes [numIndCompartmentVolumes] // 32 + * initCompartmentVolumes [numInitCompartmentVolumes] // 33 + * initFloatingSpeciesAmounts [numInitFloatingSpecies] // 34 + * boundarySpeciesAmounts [numIndBoundarySpecies] // 35 + * initBoundarySpeciesAmounts [numInitBoundarySpecies] // 36 + * globalParameters [numIndGlobalParameters] // 37 + * initGlobalParameters [numInitGlobalParameters] // 38 + * reactionRates [numReactions] // 39 * - * rateRuleValues [numRateRules] // 39 - * floatingSpeciesAmounts [numIndFloatingSpecies] // 40 + * rateRuleValues [numRateRules] // 40 + * floatingSpeciesAmounts [numIndFloatingSpecies] // 41 */ /** diff --git a/source/llvm/LLVMModelDataSymbols.cpp b/source/llvm/LLVMModelDataSymbols.cpp index daed6d5491..31e635c0bb 100644 --- a/source/llvm/LLVMModelDataSymbols.cpp +++ b/source/llvm/LLVMModelDataSymbols.cpp @@ -72,34 +72,35 @@ static const char* modelDataFieldsNames[] = { "Stoichiometry", // 13 "RandomPtr", // 14 "NumEvents", // 15 - "StateVectorSize", // 16 - "StateVector", // 17 - "StateVectorRate", // 18 - "RateRuleRates", // 19 - "FloatingSpeciesAmountRates", // 20 - - "CompartmentVolumesAlias", // 21 - "InitCompartmentVolumesAlias", // 22 - "InitFloatingSpeciesAmountsAlias", // 23 - "BoundarySpeciesAmountsAlias", // 24 - "InitBoundarySpeciesAmountsAlias", // 25 - "GlobalParametersAlias", // 26 - "InitGlobalParametersAlias", // 27 - "ReactionRatesAlias", // 28 - - "RateRuleValuesAlias", // 29 - "FloatingSpeciesAmountsAlias", // 30 - - "CompartmentVolumes", // 31 - "InitCompartmentVolumes", // 32 - "InitFloatingSpeciesAmounts", // 33 - "BoundarySpeciesAmounts", // 34 - "InitBoundarySpeciesAmounts", // 35 - "GlobalParameters", // 36 - "InitGlobalParameters", // 37 - "ReactionRates", // 38 - "NotSafe_RateRuleValues", // 39 - "NotSafe_FloatingSpeciesAmounts" // 40 + "NumPiecewiseTriggers", // 16 + "StateVectorSize", // 17 + "StateVector", // 18 + "StateVectorRate", // 19 + "RateRuleRates", // 20 + "FloatingSpeciesAmountRates", // 21 + + "CompartmentVolumesAlias", // 22 + "InitCompartmentVolumesAlias", // 23 + "InitFloatingSpeciesAmountsAlias", // 24 + "BoundarySpeciesAmountsAlias", // 25 + "InitBoundarySpeciesAmountsAlias", // 26 + "GlobalParametersAlias", // 27 + "InitGlobalParametersAlias", // 28 + "ReactionRatesAlias", // 29 + + "RateRuleValuesAlias", // 30 + "FloatingSpeciesAmountsAlias", // 31 + + "CompartmentVolumes", // 32 + "InitCompartmentVolumes", // 33 + "InitFloatingSpeciesAmounts", // 34 + "BoundarySpeciesAmounts", // 35 + "InitBoundarySpeciesAmounts", // 36 + "GlobalParameters", // 37 + "InitGlobalParameters", // 38 + "ReactionRates", // 39 + "NotSafe_RateRuleValues", // 40 + "NotSafe_FloatingSpeciesAmounts" // 41 }; @@ -333,47 +334,6 @@ int LLVMModelDataSymbols::getGlobalParameterIndex( return -1; } -/* -void LLVMModelDataSymbols::initAllocModelDataBuffers(LLVMModelData& m) const -{ - // zero out the structure - LLVMModelData::init(m); - - // set the buffer sizes - m.numIndFloatingSpecies = independentFloatingSpeciesSize; - - //mData.numDependentSpecies = ms.mNumDependentSpecies; - m.numIndGlobalParameters = independentGlobalParameterSize; - m.numReactions = reactionsMap.size(); - m.numEvents = eventAttributes.size(); - m.numRateRules = rateRules.size(); - m.numIndCompartments = independentCompartmentSize; - m.numIndBoundarySpecies = independentBoundarySpeciesSize; - - // in certain cases, the data returned by c++ new may be alligned differently than - // malloc, so just use calloc here just to be safe, plus calloc returns zero - // initialized memory. - - m.floatingSpeciesAmountsAlias = (double*)calloc(m.numIndFloatingSpecies, sizeof(double)); - m.floatingSpeciesAmountRates = 0; - m.rateRuleValuesAlias = (double*)calloc(m.numRateRules, sizeof(double)); - m.rateRuleRates = 0; - - m.reactionRatesAlias = (double*)calloc(m.numReactions, sizeof(double)); - - m.globalParametersAlias = (double*)calloc(m.numIndGlobalParameters, sizeof(double)); - m.compartmentVolumesAlias = (double*)calloc(m.numIndCompartments, sizeof(double)); - - m.boundarySpeciesAmountsAlias = (double*)calloc(m.numIndBoundarySpecies, sizeof(double)); - - - // allocate the stoichiometry matrix - m.stoichiometry = rr::csr_matrix_new(m.numIndFloatingSpecies, getReactionSize(), - stoichRowIndx, stoichColIndx, std::vector(stoichRowIndx.size(), 0)); -} -*/ - - const std::vector& LLVMModelDataSymbols::getStoichRowIndx() const { return stoichRowIndx; @@ -1665,17 +1625,42 @@ size_t LLVMModelDataSymbols::getEventBufferSize(size_t eventId) const bool LLVMModelDataSymbols::isNamedSpeciesReference(const std::string& id) const { - return namedSpeciesReferenceInfo.find(id) != namedSpeciesReferenceInfo.end(); + if (namedSpeciesReferenceInfo.find(id) != namedSpeciesReferenceInfo.end()) { + return true; + } + if (id.find(":") != string::npos) { + //LS: This whole section is for when we ask for rateOf(specid). + string rxnid = id.substr(0, id.find(":")); + string specid = id.substr(id.find(":") + 1, id.length()); + if (getReactionIndex(rxnid) != -1 && getFloatingSpeciesIndex(specid) != -1) { + return true; + } + } + return false; } const LLVMModelDataSymbols::SpeciesReferenceInfo& - LLVMModelDataSymbols::getNamedSpeciesReferenceInfo(const std::string& id) const + LLVMModelDataSymbols::getNamedSpeciesReferenceInfo(const std::string& id) { StringRefInfoMap::const_iterator i = namedSpeciesReferenceInfo.find(id); if (i != namedSpeciesReferenceInfo.end()) { return i->second; } + else if (id.find(":") != string::npos) { + //LS: This whole section is for when we ask for rateOf(specid). + string rxnid = id.substr(0, id.find(":")); + string specid = id.substr(id.find(":") + 1, id.length()); + int rxnIdx = getReactionIndex(rxnid); + int speciesIdx = getFloatingSpeciesIndex(specid); + if (rxnIdx == -1 || speciesIdx == -1) { + throw_llvm_exception(id + " is not a named SpeciesReference: '" + rxnid + "' and '" + specid + "' are not a valid combination of reaction and species."); + } + SpeciesReferenceInfo info = + { static_cast(speciesIdx), static_cast(rxnIdx), Product, rxnid }; + namedSpeciesReferenceInfo[id] = info; + return namedSpeciesReferenceInfo[id]; + } else { throw_llvm_exception(id + " is not a named SpeciesReference"); diff --git a/source/llvm/LLVMModelDataSymbols.h b/source/llvm/LLVMModelDataSymbols.h index b55a1c2258..dfef0fd66b 100644 --- a/source/llvm/LLVMModelDataSymbols.h +++ b/source/llvm/LLVMModelDataSymbols.h @@ -51,34 +51,35 @@ enum ModelDataFields { Stoichiometry, // 13 RandomPtr, // 14 NumEvents, // 15 - StateVectorSize, // 16 - StateVector, // 17 - StateVectorRate, // 18 - RateRuleRates, // 19 - FloatingSpeciesAmountRates, // 20 - - CompartmentVolumesAlias, // 21 - InitCompartmentVolumesAlias, // 22 - InitFloatingSpeciesAmountsAlias, // 23 - BoundarySpeciesAmountsAlias, // 24 - InitBoundarySpeciesAmountsAlias, // 25 - GlobalParametersAlias, // 26 - InitGlobalParametersAlias, // 27 - ReactionRatesAlias, // 28 - - RateRuleValuesAlias, // 29 - FloatingSpeciesAmountsAlias, // 30 - - CompartmentVolumes, // 31 - InitCompartmentVolumes, // 32 - InitFloatingSpeciesAmounts, // 33 - BoundarySpeciesAmounts, // 34 - InitBoundarySpeciesAmounts, // 35 - GlobalParameters, // 36 - InitGlobalParameters, // 37 - ReactionRates, // 38 - NotSafe_RateRuleValues, // 39 - NotSafe_FloatingSpeciesAmounts, // 40 + NumPiecewiseTriggers, // 16 + StateVectorSize, // 17 + StateVector, // 18 + StateVectorRate, // 19 + RateRuleRates, // 20 + FloatingSpeciesAmountRates, // 21 + + CompartmentVolumesAlias, // 22 + InitCompartmentVolumesAlias, // 23 + InitFloatingSpeciesAmountsAlias, // 24 + BoundarySpeciesAmountsAlias, // 25 + InitBoundarySpeciesAmountsAlias, // 26 + GlobalParametersAlias, // 27 + InitGlobalParametersAlias, // 28 + ReactionRatesAlias, // 29 + + RateRuleValuesAlias, // 30 + FloatingSpeciesAmountsAlias, // 31 + + CompartmentVolumes, // 32 + InitCompartmentVolumes, // 33 + InitFloatingSpeciesAmounts, // 34 + BoundarySpeciesAmounts, // 35 + InitBoundarySpeciesAmounts, // 36 + GlobalParameters, // 37 + InitGlobalParameters, // 38 + ReactionRates, // 39 + NotSafe_RateRuleValues, // 40 + NotSafe_FloatingSpeciesAmounts, // 41 }; enum EventAtributes @@ -117,7 +118,7 @@ enum EventAtributes * * Most of the indexes used in this class are indexes into ModelData * arrays, therefore we use unsigned integers -- these are less error * prone and its easier to check if they are valid i.e. only check that - * they are less than the the array size and we do not have to check that + * they are less than the array size and we do not have to check that * it is positive. * * * All symbols from the sbml are reordered such that the independent @@ -356,7 +357,7 @@ class LLVMModelDataSymbols bool isNamedSpeciesReference(const std::string& id) const; const SpeciesReferenceInfo& getNamedSpeciesReferenceInfo( - const std::string& id) const; + const std::string& id); /******* Conserved Moiety Section ********************************************/ @@ -486,12 +487,12 @@ class LLVMModelDataSymbols size_t getEventBufferSize(size_t eventId) const; /** - * the the row indices of non-zero stoichiometry values + * the row indices of non-zero stoichiometry values */ const std::vector& getStoichRowIndx() const; /** - * the the column indices of non-zero stoichiometry values + * the column indices of non-zero stoichiometry values */ const std::vector& getStoichColIndx() const; diff --git a/source/llvm/LLVMModelGenerator.cpp b/source/llvm/LLVMModelGenerator.cpp index e4585a9cd8..217ea61abf 100644 --- a/source/llvm/LLVMModelGenerator.cpp +++ b/source/llvm/LLVMModelGenerator.cpp @@ -104,6 +104,7 @@ namespace rrllvm { dst->getEventDelayPtr = src->getEventDelayPtr; dst->eventTriggerPtr = src->eventTriggerPtr; dst->eventAssignPtr = src->eventAssignPtr; + dst->getPiecewiseTriggerPtr = src->getPiecewiseTriggerPtr; dst->evalVolatileStoichPtr = src->evalVolatileStoichPtr; dst->evalConversionFactorPtr = src->evalConversionFactorPtr; } @@ -123,6 +124,7 @@ namespace rrllvm { GetEventDelayCodeGen(context).createFunction(); EventTriggerCodeGen(context).createFunction(); EventAssignCodeGen(context).createFunction(); + GetPiecewiseTriggerCodeGen(context).createFunction(); EvalVolatileStoichCodeGen(context).createFunction(); EvalConversionFactorCodeGen(context).createFunction(); @@ -252,7 +254,7 @@ namespace rrllvm { modelGeneratorContext->getJitNonOwning()->addModule(); LLVMModelData* modelData = createModelData(modelGeneratorContext->getModelDataSymbols(), - modelGeneratorContext->getRandom()); + modelGeneratorContext->getRandom(), modelGeneratorContext->getNumPiecewiseTriggers()); uint llvmsize = ModelDataIRBuilder::getModelDataSize( modelGeneratorContext->getJitNonOwning()->getModuleNonOwning(), @@ -520,7 +522,7 @@ namespace rrllvm { if (sp) { rrLog(Logger::LOG_DEBUG) << "found a cached model for \"" << sbmlMD5 << "\""; - return new LLVMExecutableModel(sp, createModelData(*sp->symbols, sp->random)); + return new LLVMExecutableModel(sp, createModelData(*sp->symbols, sp->random, 0)); } else { rrLog(Logger::LOG_DEBUG) << "no cached model found for " << sbmlMD5 @@ -611,7 +613,7 @@ namespace rrllvm { return str; } - LLVMModelData* createModelData(const rrllvm::LLVMModelDataSymbols& symbols, const Random* random) { + LLVMModelData* createModelData(const rrllvm::LLVMModelDataSymbols& symbols, const Random* random, uint numPiecewiseTriggers) { uint modelDataBaseSize = sizeof(LLVMModelData); uint numIndCompartments = static_cast(symbols.getIndependentCompartmentSize()); @@ -659,6 +661,7 @@ namespace rrllvm { modelData->numRateRules = numRateRules; modelData->numReactions = numReactions; modelData->numEvents = static_cast(symbols.getEventAttributes().size()); + modelData->numPiecewiseTriggers = numPiecewiseTriggers; // set the aliases to the offsets uint offset = 0; diff --git a/source/llvm/LoadSymbolResolverBase.cpp b/source/llvm/LoadSymbolResolverBase.cpp index b881449e78..3746c93fe4 100644 --- a/source/llvm/LoadSymbolResolverBase.cpp +++ b/source/llvm/LoadSymbolResolverBase.cpp @@ -82,6 +82,11 @@ void LoadSymbolResolverBase::recursiveSymbolPop() symbolStack.pop_back(); } +bool LoadSymbolResolverBase::isLocalParameter(const std::string& symbol) +{ + return false; +} + void LoadSymbolResolverBase::flushCache() { symbolCache.clear(); diff --git a/source/llvm/LoadSymbolResolverBase.h b/source/llvm/LoadSymbolResolverBase.h index 4f43f78939..678a28952d 100644 --- a/source/llvm/LoadSymbolResolverBase.h +++ b/source/llvm/LoadSymbolResolverBase.h @@ -43,6 +43,8 @@ class LoadSymbolResolverBase: public LoadSymbolResolver void recursiveSymbolPop() override; + virtual bool isLocalParameter(const std::string& symbol) override; + /** * Flush the symbol cache. This is required in branches and switch blocks as * a symbol used in a previous block can not be re-used in the current block. diff --git a/source/llvm/ModelDataIRBuilder.cpp b/source/llvm/ModelDataIRBuilder.cpp index 3a18895a3b..f7ce1363e7 100644 --- a/source/llvm/ModelDataIRBuilder.cpp +++ b/source/llvm/ModelDataIRBuilder.cpp @@ -52,23 +52,22 @@ static bool isAliasOrPointer(ModelDataFields f) { /* Model Data alias are between one of the following values */ /* - StateVector, // 13 - StateVectorRate, // 14 - RateRuleRates, // 15 - FloatingSpeciesAmountRates, // 16 - - CompartmentVolumesAlias, // 17 - CompartmentVolumesInitAlias, // 18 - FloatingSpeciesAmountsInitAlias, // 19 - ConservedSpeciesAmountsInitAlias, // 20 - BoundarySpeciesAmountsAlias, // 21 - BoundarySpeciesAmountsInitAlias, // 22 - GlobalParametersAlias, // 23 - GlobalParametersInitAlias, // 24 - ReactionRatesAlias, // 25 - - RateRuleValuesAlias, // 26 - FloatingSpeciesAmountsAlias, // 27 + "StateVector", // 18 + "StateVectorRate", // 19 + "RateRuleRates", // 20 + "FloatingSpeciesAmountRates", // 21 + + "CompartmentVolumesAlias", // 22 + "InitCompartmentVolumesAlias", // 23 + "InitFloatingSpeciesAmountsAlias", // 24 + "BoundarySpeciesAmountsAlias", // 25 + "InitBoundarySpeciesAmountsAlias", // 26 + "GlobalParametersAlias", // 27 + "InitGlobalParametersAlias", // 28 + "ReactionRatesAlias", // 29 + + "RateRuleValuesAlias", // 30 + "FloatingSpeciesAmountsAlias", // 31 */ return f >= StateVector && f <= FloatingSpeciesAmountsAlias; } @@ -77,17 +76,16 @@ static bool isArray(ModelDataFields f) { /* arrays are the last elements in the model data struct */ /* - CompartmentVolumes, // 28 - CompartmentVolumesInit, // 29 - FloatingSpeciesAmounts, // 30 - FloatingSpeciesAmountsInit, // 31 - ConservedSpeciesAmountsInit, // 32 - BoundarySpeciesAmounts, // 33 - BoundarySpeciesAmountsInit, // 34 - GlobalParameters, // 35 - GlobalParametersInit, // 36 - RateRuleValues, // 37 - ReactionRates, // 38 + "CompartmentVolumes", // 32 + "InitCompartmentVolumes", // 33 + "InitFloatingSpeciesAmounts", // 34 + "BoundarySpeciesAmounts", // 35 + "InitBoundarySpeciesAmounts", // 36 + "GlobalParameters", // 37 + "InitGlobalParameters", // 38 + "ReactionRates", // 39 + "NotSafe_RateRuleValues", // 40 + "NotSafe_FloatingSpeciesAmounts" // 41 */ return f >= CompartmentVolumes && f <= NotSafe_FloatingSpeciesAmounts; @@ -687,34 +685,35 @@ llvm::StructType *ModelDataIRBuilder::createModelDataStructType(llvm::Module *mo elements.push_back(csrSparsePtrType); // 13 dcsr_matrix stoichiometry; elements.push_back(voidPtrType); // 14 void* random; elements.push_back(int32Type); // 15 int numEvents; - elements.push_back(int32Type); // 16 int stateVectorSize; - elements.push_back(doublePtrType); // 17 double* stateVector; - elements.push_back(doublePtrType); // 18 double* stateVectorRate; - elements.push_back(doublePtrType); // 19 double* rateRuleRates; - elements.push_back(doublePtrType); // 20 double* floatingSpeciesAmountRates; - - elements.push_back(doublePtrType); // 21 double* compartmentVolumesAlias; - elements.push_back(doublePtrType); // 22 double* compartmentVolumesInitAlias; - elements.push_back(doublePtrType); // 23 double* floatingSpeciesAmountsInitAlias - elements.push_back(doublePtrType); // 24 double* boundarySpeciesAmountsAlias; - elements.push_back(doublePtrType); // 25 double* boundarySpeciesAmountsInitAlias; - elements.push_back(doublePtrType); // 26 double* globalParametersAlias - elements.push_back(doublePtrType); // 27 double* globalParametersInitAlias - elements.push_back(doublePtrType); // 28 double* reactionRatesAlias - - elements.push_back(doublePtrType); // 29 double* rateRuleValuesAlias - elements.push_back(doublePtrType); // 30 double* floatingSpeciesAmountsAlias - - elements.push_back(ArrayType::get(doubleType, numIndCompartments)); // 31 CompartmentVolumes - elements.push_back(ArrayType::get(doubleType, numInitCompartments)); // 32 initCompartmentVolumes - elements.push_back(ArrayType::get(doubleType, numInitFloatingSpecies)); // 33 initFloatingSpeciesAmounts - elements.push_back(ArrayType::get(doubleType, numIndBoundarySpecies)); // 34 boundarySpeciesAmounts - elements.push_back(ArrayType::get(doubleType, numInitBoundarySpecies)); // 35 initBoundarySpeciesAmounts - elements.push_back(ArrayType::get(doubleType, numIndGlobalParameters)); // 36 globalParameters - elements.push_back(ArrayType::get(doubleType, numInitGlobalParameters));// 37 initGlobalParameters - elements.push_back(ArrayType::get(doubleType, numReactions)); // 38 reactionRates - elements.push_back(ArrayType::get(doubleType, numRateRules)); // 39 rateRuleValues - elements.push_back(ArrayType::get(doubleType, numIndFloatingSpecies)); // 40 floatingSpeciesAmounts + elements.push_back(int32Type); // 16 int numPiecewiseTriggers; + elements.push_back(int32Type); // 17 int stateVectorSize; + elements.push_back(doublePtrType); // 18 double* stateVector; + elements.push_back(doublePtrType); // 19 double* stateVectorRate; + elements.push_back(doublePtrType); // 20 double* rateRuleRates; + elements.push_back(doublePtrType); // 21 double* floatingSpeciesAmountRates; + + elements.push_back(doublePtrType); // 22 double* compartmentVolumesAlias; + elements.push_back(doublePtrType); // 23 double* compartmentVolumesInitAlias; + elements.push_back(doublePtrType); // 24 double* floatingSpeciesAmountsInitAlias + elements.push_back(doublePtrType); // 25 double* boundarySpeciesAmountsAlias; + elements.push_back(doublePtrType); // 26 double* boundarySpeciesAmountsInitAlias; + elements.push_back(doublePtrType); // 27 double* globalParametersAlias + elements.push_back(doublePtrType); // 28 double* globalParametersInitAlias + elements.push_back(doublePtrType); // 29 double* reactionRatesAlias + + elements.push_back(doublePtrType); // 30 double* rateRuleValuesAlias + elements.push_back(doublePtrType); // 31 double* floatingSpeciesAmountsAlias + + elements.push_back(ArrayType::get(doubleType, numIndCompartments)); // 32 CompartmentVolumes + elements.push_back(ArrayType::get(doubleType, numInitCompartments)); // 33 initCompartmentVolumes + elements.push_back(ArrayType::get(doubleType, numInitFloatingSpecies)); // 34 initFloatingSpeciesAmounts + elements.push_back(ArrayType::get(doubleType, numIndBoundarySpecies)); // 35 boundarySpeciesAmounts + elements.push_back(ArrayType::get(doubleType, numInitBoundarySpecies)); // 36 initBoundarySpeciesAmounts + elements.push_back(ArrayType::get(doubleType, numIndGlobalParameters)); // 37 globalParameters + elements.push_back(ArrayType::get(doubleType, numInitGlobalParameters));// 38 initGlobalParameters + elements.push_back(ArrayType::get(doubleType, numReactions)); // 39 reactionRates + elements.push_back(ArrayType::get(doubleType, numRateRules)); // 40 rateRuleValues + elements.push_back(ArrayType::get(doubleType, numIndFloatingSpecies)); // 41 floatingSpeciesAmounts // creates a named struct, // the act of creating a named struct should @@ -956,377 +955,6 @@ llvm::Function* LLVMModelDataIRBuilderTesting::getDispIntDecl(llvm::Module* modu } - - -void LLVMModelDataIRBuilderTesting::test(Module *module, IRBuilder<> *build, - ExecutionEngine * engine) -{ -// TheModule = module; -// context = &module->getContext(); -// rr::builder = build; -// TheExecutionEngine = engine; -// createDispPrototypes(); -// -// callDispInt(23); -// -// structSetProto(); -// -// TestStruct s = -// { 0 }; -// -// s.var3 = (char*) calloc(50, sizeof(char)); -// s.var5 = (char*) calloc(50, sizeof(char)); -// -// sprintf(s.var5, "1234567890"); -// -// char* p = (char*) calloc(50, sizeof(char)); -// -// dispStruct(&s); -// -// callStructSet(&s, 314, 3.14, 2.78, p, 0, 0); -// -// printf("p: %s\n", p); -// -// dispStruct(&s); -// -// printf("new var5: "); -// for (int i = 0; i < 10; i++) -// { -// printf("{i:%i,c:%i}, ", i, s.var5[i]); -// } -// printf("\n"); -} - - - - -// -// -//static Module *TheModule; -//static IRBuilder<> *builder; -//static ExecutionEngine *TheExecutionEngine; -//static LLVMContext *context; -// -//static LLVMContext &getContext() { -// return *context; -//} -// -//Value *storeInStruct(IRBuilder<> *B, Value *S, Value *V, unsigned index, const char* s = ""); -//Value *loadFromStruct(IRBuilder<> *B, Value *S, unsigned index, const char* s = ""); -// -//struct TestStruct { -// int var0; -// double var1; -// double var2; -// char* var3; -// double var4; -// char* var5; -//}; -// -//StructType* getTestStructType() { -// // static StructType * create (ArrayRef< Type * > Elements, StringRef Name, bool isPacked=false) -// std::vector elements; -// elements.push_back(Type::getInt32Ty(getContext())); // int var0; -// elements.push_back(Type::getDoubleTy(getContext())); // double var1; -// elements.push_back(Type::getDoubleTy(getContext())); // double var2; -// //elements.push_back(ArrayType::get(Type::getInt8Ty(getContext()), 0)); // char* var3; -// elements.push_back(Type::getInt8PtrTy(getContext())); // char* var5; -// elements.push_back(Type::getDoubleTy(getContext())); // double var4; -// //elements.push_back(ArrayType::get(Type::getInt8Ty(getContext()), 0)); // char[] var5; -// elements.push_back(Type::getInt8PtrTy(getContext())); // char* var5; // char* var5; -// -// -// -// -// StructType *s = StructType::create(elements, "TestStruct"); -// -// const DataLayout &dl = *TheExecutionEngine->getDataLayout(); -// -// size_t llvm_size = dl.getTypeStoreSize(s); -// -// printf("TestStruct size: %i, , LLVM Size: %i\n", sizeof(TestStruct), llvm_size); -// -// -// -// return s; -//} -// -//Value *getVar5(Value *testStructPtr, int index) { -// GetElementPtrInst *gep = dyn_cast(builder->CreateStructGEP(testStructPtr, 5, "var5_gep")); -// LoadInst * var5_load = builder.CreateLoad(gep, "var5_load"); -// //Value *var5_load = loadFromStruct(Builder, testStructPtr, 5, "var5"); -// gep = dyn_cast(builder->CreateConstGEP1_32(var5_load, index, "var5_elem_gep")); -// return builder.CreateLoad(gep, "var5_elem"); -//} -// -// -//static void dispStruct(TestStruct *s) -//{ -// printf("%s {\n", __PRETTY_FUNCTION__); -// printf("\tvar0: %i\n", s->var0); -// printf("\tvar1: %f\n", s->var1); -// printf("\tvar2: %f\n", s->var2); -// printf("\tvar3: %s\n", s->var3); -// printf("\tvar4: %f\n", s->var4); -// printf("\tvar5: %s\n", s->var5); -// printf("}\n"); -// -// -//} -// -//static Function* dispStructProto() { -// -// StructType *structType = getTestStructType(); -// PointerType *structTypePtr = llvm::PointerType::get(structType, 0); -// std::vector args(1, structTypePtr); -// FunctionType *funcType = FunctionType::get(Type::getVoidTy(getContext()), args, false); -// Function *f = Function::Create(funcType, Function::InternalLinkage, "dispStruct", TheModule); -// return f; -//} -// -//static void dispInt(int i) { -// printf("%s, %i\n", __PRETTY_FUNCTION__, i); -//} -// - -// -//static void callDispInt(int val) { -// // Evaluate a top-level expression into an anonymous function. -// -// Function *func = TheExecutionEngine->FindFunctionNamed("dispInt"); -// -// // JIT the function, returning a function pointer. -// void *fptr = TheExecutionEngine->getPointerToFunction(func); -// -// // Cast it to the right type (takes no arguments, returns a double) so we -// // can call it as a native function. -// void (*fp)(int) = (void (*)(int))fptr; -// -// fp(val); -//} -// -// -// - -// -//static void dispCharStar(char* p) { -// printf("%s: %s\n", __PRETTY_FUNCTION__, p); -//} -// -//static Function* dispCharStarProto() { -// Function *f = TheModule->getFunction("dispCharStar"); -// if (f == 0) { -// std::vectorargs(1, Type::getInt8PtrTy(getContext())); -// FunctionType *funcType = FunctionType::get(Type::getVoidTy(getContext()), args, false); -// f = Function::Create(funcType, Function::InternalLinkage, "dispCharStar", TheModule); -// } -// return f; -//} -// -//static void dispChar(char p) { -// printf("%s as char: %c\n", __PRETTY_FUNCTION__, p); -// printf("%s as int: %i\n", __PRETTY_FUNCTION__, (int)p); -//} -// -//static Function* dispCharProto() { -// Function *f = TheModule->getFunction("dispChar"); -// if (f == 0) { -// std::vectorargs(1, Type::getInt8Ty(getContext())); -// FunctionType *funcType = FunctionType::get(Type::getVoidTy(getContext()), args, false); -// f = Function::Create(funcType, Function::InternalLinkage, "dispChar", TheModule); -// } -// return f; -//} -// -//static void createDispPrototypes() { -// TheExecutionEngine->addGlobalMapping(dispStructProto(), (void*)dispStruct); -// TheExecutionEngine->addGlobalMapping(dispIntProto(), (void*)dispInt); -// TheExecutionEngine->addGlobalMapping(dispDoubleProto(), (void*)dispDouble); -// TheExecutionEngine->addGlobalMapping(dispCharStarProto(), (void*)dispCharStar); -// TheExecutionEngine->addGlobalMapping(dispCharProto(), (void*)dispChar); -//} -// -//static void structSetBody(Function *func) { -// // Create a new basic block to start insertion into. -// BasicBlock *BB = BasicBlock::Create(getContext(), "entry", func); -// builder.SetInsertPoint(BB); -// -// LLVMContext &context = getContext(); -// -// std::vector args; -// -// // Set names for all arguments. -// unsigned idx = 0; -// for (Function::arg_iterator ai = func->arg_begin(); ai != func->arg_end(); -// ++ai, ++idx) { -// -// args.push_back(ai); -// //ai->setName(Args[Idx]); -// -// // Add arguments to variable symbol table. -// //NamedValues[Args[Idx]] = AI; -// } -// -// Value *arg = ConstantInt::get(Type::getInt32Ty(getContext()), 16); -// -// builder.CreateCall(dispIntProto(), args[1], ""); -// -// builder.CreateCall(dispDoubleProto(), args[2], ""); -// -// //builder->CreateCall(dispDoubleProto(), args[3], ""); -// -// builder.CreateCall(dispCharStarProto(), args[4], ""); -// -// storeInStruct(builder, args[0], args[1], 0); -// storeInStruct(builder, args[0], args[2], 1); -// storeInStruct(builder, args[0], args[3], 2); -// storeInStruct(builder, args[0], args[5], 4); -// -// -// -// Value *var3 = loadFromStruct(builder, args[0], 3, "var3"); -// -// //Value *var3 = args[4]; -// -// Type *t = var3->getType(); -// -// var3->dump(); -// t->dump(); -// -// -// Value *c = ConstantInt::get(Type::getInt8Ty(getContext()), 's'); -// Value *gep = builder.CreateConstGEP1_32(var3, 0, "array gep"); -// builder.CreateStore(c, gep); -// gep = builder.CreateConstGEP1_32(var3, 1, "array gep"); -// builder.CreateStore(c, gep); -// gep = builder.CreateConstGEP1_32(var3, 2, "array gep"); -// builder.CreateStore(c, gep); -// -// -// c = ConstantInt::get(Type::getInt8Ty(getContext()), '2'); -// -// Value *var5_2 = getVar5(args[0], 2); -// -// var5_2->dump(); -// -// builder.CreateCall(dispCharProto(), var5_2); -// -// -// -// builder.CreateCall(dispCharStarProto(), args[4], ""); -// -// -// -// -// // Finish off the function. -// builder.CreateRetVoid(); -// -// // Validate the generated code, checking for consistency. -// verifyFunction(*func); -//} -// -//static Function* structSetProto() { -// Function *f = TheModule->getFunction("structSet"); -// -// if (f == 0) { -// std::vector args; -// StructType *structType = getTestStructType(); -// PointerType *structTypePtr = llvm::PointerType::get(structType, 0); -// args.push_back(structTypePtr); -// args.push_back(Type::getInt32Ty(getContext())); // int var0; -// args.push_back(Type::getDoubleTy(getContext())); // double var1; -// args.push_back(Type::getDoubleTy(getContext())); // double var2; -// args.push_back(Type::getInt8PtrTy(getContext())); // char* var3; -// args.push_back(Type::getDoubleTy(getContext())); // double var4; -// args.push_back(Type::getInt8PtrTy(getContext())); // char* var5; -// -// FunctionType *funcType = FunctionType::get(Type::getVoidTy(getContext()), args, false); -// f = Function::Create(funcType, Function::InternalLinkage, "structSet", TheModule); -// -// -// structSetBody(f); -// -// } -// return f; -//} -// -//static void callStructSet(TestStruct *s, int var0, double var1, double var2, -// char* var3, double var4, char* var5) -//{ -// Function *func = TheExecutionEngine->FindFunctionNamed("structSet"); -// -// // JIT the function, returning a function pointer. -// void *fptr = TheExecutionEngine->getPointerToFunction(func); -// -// // Cast it to the right type (takes no arguments, returns a double) so we -// // can call it as a native function. -// void (*fp)(TestStruct*, int, double, double, char*, double, char*) = -// (void (*)(TestStruct*, int, double, double,char*, double, char*))fptr; -// -// fp(s, var0, var1, var2, var3, var4, var5); -// -// -// -//} -// -// -// -//// Store V in structure S element index -//Value *storeInStruct(IRBuilder<> *B, Value *S, Value *V, unsigned index, const char* s) -//{ -// return B->CreateStore(V, B->CreateStructGEP(S, index, s), s); -//} -// -//// Store V in structure S element index -//Value *loadFromStruct(IRBuilder<> *B, Value *S, unsigned index, const char* s) -//{ -// Value *gep = B->CreateStructGEP(S, index, s); -// return B->CreateLoad(gep, s); -//} -// -//void dispIntrinsics(Intrinsic::ID id) { -// Intrinsic::IITDescriptor ids[8]; -// SmallVector table; -// Intrinsic::getIntrinsicInfoTableEntries(id, table); -// -// printf("table: %i\n", table.size()); -// -// for(unsigned i = 0; i < table.size(); i++) { -// ids[i] = table[i]; -// } -// -// for(unsigned i = 0; i < table.size(); i++) { -// Intrinsic::IITDescriptor::IITDescriptorKind kind = ids[i].Kind; -// Intrinsic::IITDescriptor::ArgKind argKind = ids[i].getArgumentKind(); -// unsigned argNum = ids[i].getArgumentNumber(); -// -// printf("kind: %i, argKind: %i, argNum: %i\n", kind, argKind, argNum); -// -// } -// -// printf("foo\n"); -// -// std::vector args(table.size() - 1, Type::getDoubleTy(getContext())); -// -// printf("intrinsic name: %s\n", getName(id, args).c_str()); -// -// Function *decl = Intrinsic::getDeclaration(TheModule, id, args); -// decl->dump(); -// -// FunctionType *func = Intrinsic::getType(getContext(), id, args); -// -// func->dump(); -// -// printf("done\n"); -// -// -// -// -// -// //func->dump(); -//} - - } /* namespace rr */ diff --git a/source/llvm/ModelDataIRBuilder.h b/source/llvm/ModelDataIRBuilder.h index fd9539a71d..c8d7345a4c 100644 --- a/source/llvm/ModelDataIRBuilder.h +++ b/source/llvm/ModelDataIRBuilder.h @@ -380,9 +380,6 @@ namespace rrllvm std::pair createFloatingSpeciesAccessors( llvm::Module* module, const std::string id); - void test(llvm::Module* module, llvm::IRBuilder<>* build, - llvm::ExecutionEngine* engine); - static llvm::Function* getDispIntDecl(llvm::Module* module); llvm::CallInst* createDispInt(llvm::Value* intVal); diff --git a/source/llvm/ModelDataSymbolResolver.cpp b/source/llvm/ModelDataSymbolResolver.cpp index fe6677db5c..fcf957e45d 100644 --- a/source/llvm/ModelDataSymbolResolver.cpp +++ b/source/llvm/ModelDataSymbolResolver.cpp @@ -148,10 +148,12 @@ namespace rrllvm return cacheValue(symbol, args, mdbuilder.createRateRuleValueLoad(symbol)); } + LLVMModelDataSymbols* modelDataSymbolsPtr = const_cast(&modelDataSymbols); + if (modelDataSymbols.isNamedSpeciesReference(symbol)) { const LLVMModelDataSymbols::SpeciesReferenceInfo& info = - modelDataSymbols.getNamedSpeciesReferenceInfo(symbol); + modelDataSymbolsPtr->getNamedSpeciesReferenceInfo(symbol); Value* value = mdbuilder.createStoichiometryLoad(info.row, info.column, symbol); @@ -265,6 +267,8 @@ namespace rrllvm // at this point, we have already taken care of the species amount / // conc conversion, rest are just plain stores. + LLVMModelDataSymbols* modelDataSymbolsPtr = const_cast(&modelDataSymbols); + if (modelDataSymbols.hasRateRule(symbol)) { return mdbuilder.createRateRuleValueStore(symbol, value); @@ -283,7 +287,7 @@ namespace rrllvm else if (modelDataSymbols.isNamedSpeciesReference(symbol)) { const LLVMModelDataSymbols::SpeciesReferenceInfo& info = - modelDataSymbols.getNamedSpeciesReferenceInfo(symbol); + modelDataSymbolsPtr->getNamedSpeciesReferenceInfo(symbol); if (info.type == LLVMModelDataSymbols::MultiReactantProduct) { diff --git a/source/llvm/ModelGeneratorContext.cpp b/source/llvm/ModelGeneratorContext.cpp index e19cd62161..f16edaf5bc 100644 --- a/source/llvm/ModelGeneratorContext.cpp +++ b/source/llvm/ModelGeneratorContext.cpp @@ -28,9 +28,11 @@ #include "llvm/IR/DataLayout.h" #include "llvm/Support/TargetRegistry.h" #include "llvm/Support/DynamicLibrary.h" -//#include "llvm/ExecutionEngine/OrcMCJITReplacement.h" #include "llvm/ExecutionEngine/SectionMemoryManager.h" +#include +#include + #include "llvm/ExecutionEngine/ExecutionEngine.h" @@ -108,6 +110,7 @@ namespace rrllvm { : ownedDoc(nullptr), doc(nullptr), + mPiecewiseTriggers(), symbols(nullptr), options(options), random(nullptr), @@ -161,6 +164,8 @@ namespace rrllvm { model = doc->getModel(); + addAllPiecewiseTriggers(model); + symbols = new LLVMModelDataSymbols(doc->getModel(), options); modelSymbols = std::make_unique(getModel(), *symbols); @@ -207,60 +212,155 @@ namespace rrllvm { ModelGeneratorContext::ModelGeneratorContext() : ownedDoc(createEmptyDocument()), doc(ownedDoc), + mPiecewiseTriggers(), symbols(new LLVMModelDataSymbols(doc->getModel(), 0)), modelSymbols(new LLVMModelSymbols(getModel(), *symbols)), options(0) -// functionPassManager(0) { // initialize LLVM // TODO check result InitializeNativeTarget(); InitializeNativeTargetAsmPrinter(); InitializeNativeTargetAsmParser(); - -// context = std::unique_ptr(); -// // Make the module, which holds all the code. -// module_owner = std::unique_ptr(new Module("LLVM Module", *context)); -// module = module_owner.get(); -// -// builder = std::make_unique>(*context); -// -// errString = new std::string(); -// -// addGlobalMappings(); -// -// EngineBuilder engineBuilder(std::move(module_owner)); -// //engineBuilder.setEngineKind(EngineKind::JIT); -// engineBuilder.setErrorStr(errString); -// executionEngine = std::unique_ptr(engineBuilder.create()); - } Random *ModelGeneratorContext::getRandom() const { return random; } + void ModelGeneratorContext::addAllPiecewiseTriggers(const libsbml::Model* model) + { + clearPiecewiseTriggers(); + + //If there are piecewise functions in function definitions, we need to convert the document to not have function definitions any more (sigh) + for (size_t fd = 0; fd < model->getNumFunctionDefinitions(); fd++) { + if (containsPiecewise(model->getFunctionDefinition(fd)->getMath())) { + libsbml::SBMLFunctionDefinitionConverter converter; + SBMLDocument doc(model->getLevel(), model->getVersion()); + doc.setModel(model); + converter.setDocument(&doc); + if (converter.convert() == LIBSBML_OPERATION_SUCCESS) { + return addAllPiecewiseTriggers(doc.getModel()); + } + //Otherwise, I suppose we just continue? + rrLog(Logger::LOG_WARNING) << "A piecewise function was discovered in a function definition, but we were unable to convert the document to remove function definitions. Any transitions in those piecewise functions may not be noticed by the simulator."; + } + } + + //Only worry about piecewise functions in continuous contexts: initial assignments and event assignments are fine, because we don't need to stop simulation to check things; the assignment is already happening when simulation is stopped. + for (size_t r = 0; r < model->getNumRules(); r++) { + //We don't care whether the rule is assignment, rate, or even algebraic; we're just looking for piecewise functions. + const libsbml::Rule* rule = model->getRule(r); + const ASTNode* node = rule->getMath(); + addPiecewiseTriggersFrom(node); + } + for (size_t r = 0; r < model->getNumReactions(); r++) { + const libsbml::Reaction* rxn = model->getReaction(r); + if (rxn->isSetKineticLaw()) { + addPiecewiseTriggersFrom(rxn->getKineticLaw()->getMath()); + } + } + //Event triggers themselves might have piecewise arguments in them: + for (size_t e = 0; e < model->getNumEvents(); e++) { + const libsbml::Event* event = model->getEvent(e); + if (event->isSetTrigger()) { + addPiecewiseTriggersFrom(event->getTrigger()->getMath()); + } + } + } + + void ModelGeneratorContext::addPiecewiseTriggersFrom(const libsbml::ASTNode* node) + { + if (node == NULL) { + return; + } + if (node->getType() == libsbml::AST_FUNCTION_PIECEWISE) { + // The format of the piecewise function is (val1, cond1, val2, cond2, ..., val_otherwise). They're checked sequentially, so if cond1 is true, val1 is always returned, and if cond1 is false but cond2 is true, val2 is returned, etc. This means the value returned is dependent on the following conditions: + // cond1 + // cond2 && (!cond1) + // cond3 && (!cond1 && !cond2) + // ... + // !cond1 && !cond2 && !cond3 && ... + + // None of the values need to be checked, just the conditions. Hence: + std::vector conditions; + for (size_t c = 1; c < node->getNumChildren(); c += 2) { + conditions.push_back(node->getChild(c)); + } + for (size_t c = 0; c < conditions.size(); c++) { + std::vector nots; + for (size_t n = 0; n < c; n++) { + nots.push_back(conditions[n]); + } + if (nots.size() > 0) { + libsbml::ASTNode* root = new libsbml::ASTNode(libsbml::AST_LOGICAL_AND); + root->addChild(conditions[c]->deepCopy()); + for (size_t n = 0; n < nots.size(); n++) { + libsbml::ASTNode* child_not = new libsbml::ASTNode(libsbml::AST_LOGICAL_NOT); + child_not->addChild(nots[n]->deepCopy()); + root->addChild(child_not); + } + mPiecewiseTriggers.push_back(root); + } + else { + mPiecewiseTriggers.push_back(conditions[c]->deepCopy()); + } + } + //If there's an 'otherwise', we need an all-nots trigger: + if (node->getNumChildren() % 2 == 1) { + libsbml::ASTNode* root = new libsbml::ASTNode(libsbml::AST_LOGICAL_AND); + for (size_t c = 0; c < conditions.size(); c++) { + libsbml::ASTNode* child_not = new libsbml::ASTNode(libsbml::AST_LOGICAL_NOT); + child_not->addChild(conditions[c]->deepCopy()); + root->addChild(child_not); + } + mPiecewiseTriggers.push_back(root); + } + } + + //Not 'else', since a piecewise function can have piecewise children. + for (size_t c = 0; c < node->getNumChildren(); c++) { + const ASTNode* child = node->getChild(c); + addPiecewiseTriggersFrom(child); + } + } + + bool ModelGeneratorContext::containsPiecewise(const libsbml::ASTNode* node) + { + if (node == NULL) { + return false; + } + if (node->getType() == libsbml::AST_FUNCTION_PIECEWISE) { + return true; + } + + for (size_t c = 0; c < node->getNumChildren(); c++) { + if (containsPiecewise(node->getChild(c))) { + return true; + } + } + + return false; + } + + void ModelGeneratorContext::clearPiecewiseTriggers() + { + for (size_t pt = 0; pt < mPiecewiseTriggers.size(); pt++) { + delete mPiecewiseTriggers[pt]; + } + mPiecewiseTriggers.clear(); + } + void ModelGeneratorContext::cleanup() { delete ownedDoc; ownedDoc = 0; -// delete errString; errString = 0; + clearPiecewiseTriggers(); } - ModelGeneratorContext::~ModelGeneratorContext() { cleanup(); } -//llvm::LLVMContext &ModelGeneratorContext::getContext() const -//{ -// return *context; -//} -// -//llvm::ExecutionEngine &ModelGeneratorContext::getExecutionEngine() const -//{ -// return *executionEngine; -//} - const LLVMModelDataSymbols &ModelGeneratorContext::getModelDataSymbols() const { return *symbols; } @@ -269,26 +369,23 @@ namespace rrllvm { return doc; } - Jit *ModelGeneratorContext::getJitNonOwning() const { return jit.get(); } - const libsbml::Model *ModelGeneratorContext::getModel() const { return doc->getModel(); } + const std::vector* ModelGeneratorContext::getPiecewiseTriggers() const + { + return &mPiecewiseTriggers; + } -//llvm::Module *ModelGeneratorContext::getModule() const -//{ -// return jit->getModuleNonOwning(); -//} -// -//llvm::IRBuilder<> &ModelGeneratorContext::getBuilder() const -//{ -// return *jit->getBuilderNonOwning(); -//} + size_t ModelGeneratorContext::getNumPiecewiseTriggers() const + { + return mPiecewiseTriggers.size(); + } void ModelGeneratorContext::transferObjectsToResources(std::shared_ptr modelResources) { modelResources->symbols = symbols; @@ -299,14 +396,6 @@ namespace rrllvm { modelResources->random = random; random = nullptr; - -// modelResources->context = std::move(context); -// context = nullptr; -// modelResources->executionEngine = std::move(executionEngine); -// executionEngine = nullptr; -// -// modelResources->errStr = errString; -// errString = nullptr; } const LLVMModelSymbols &ModelGeneratorContext::getModelSymbols() const { @@ -321,7 +410,6 @@ namespace rrllvm { return (options & LoadSBMLOptions::LLVM_SYMBOL_CACHE) != 0; } - static SBMLDocument *checkedReadSBMLFromString(const char *xml) { SBMLDocument *doc = readSBMLFromString(xml); diff --git a/source/llvm/ModelGeneratorContext.h b/source/llvm/ModelGeneratorContext.h index b042857636..7e26483351 100644 --- a/source/llvm/ModelGeneratorContext.h +++ b/source/llvm/ModelGeneratorContext.h @@ -120,25 +120,11 @@ namespace rrllvm { const libsbml::Model *getModel() const; - Jit *getJitNonOwning() const; + const std::vector* getPiecewiseTriggers() const; - /** - * nearly all llvm functions expect a pointer to module, so we define this - * as a pointer type instead of reference, even though we gaurantee this to - * be non-null - */ -// llvm::Module *getModule() const; + size_t getNumPiecewiseTriggers() const; - /** - * if optimization is enabled, this gets the function pass - * manager loaded with all the requested optimizers. - * NULL if no optimization is specified. - * - * @return a borrowed reference that is owned by ModelGeneratorContext - */ -// llvm::legacy::FunctionPassManager* getFunctionPassManager() const; - -// llvm::IRBuilder<> &getBuilder() const; + Jit *getJitNonOwning() const; /** * A lot can go wrong in the process of generating a model from an sbml doc. @@ -175,7 +161,7 @@ namespace rrllvm { /** * these point to the same location, ownedDoc is set if we create the doc, - * otherwise its nullptr, meaning we're borrowing the the doc. + * otherwise its nullptr, meaning we're borrowing the doc. */ libsbml::SBMLDocument *ownedDoc; @@ -184,6 +170,8 @@ namespace rrllvm { */ const libsbml::SBMLDocument *doc; + std::vector mPiecewiseTriggers; + const LLVMModelDataSymbols *symbols; /** @@ -232,6 +220,26 @@ namespace rrllvm { */ std::unique_ptr moietyConverter; + /** + * Get all transitions from any piecewise equations in the model. + */ + void addAllPiecewiseTriggers(const libsbml::Model* model); + + /** + * Get all transitions from any piecewise equations in the model. + */ + void addPiecewiseTriggersFrom(const libsbml::ASTNode* node); + + /** + * Get all transitions from any piecewise equations in the model. + */ + bool containsPiecewise(const libsbml::ASTNode* node); + + /** + * Delete any piecewise trigger nodes (which we own) + */ + void clearPiecewiseTriggers(); + /** * free any memory this class allocated. */ @@ -239,7 +247,7 @@ namespace rrllvm { }; - LLVMModelData *createModelData(const rrllvm::LLVMModelDataSymbols &symbols, const Random *random); + LLVMModelData *createModelData(const rrllvm::LLVMModelDataSymbols &symbols, const Random *random, uint numPiecewiseTriggers); } /* namespace rr */ diff --git a/source/llvm/ModelInitialValueSymbolResolver.cpp b/source/llvm/ModelInitialValueSymbolResolver.cpp index 6241b4a858..351f50aa88 100644 --- a/source/llvm/ModelInitialValueSymbolResolver.cpp +++ b/source/llvm/ModelInitialValueSymbolResolver.cpp @@ -219,6 +219,36 @@ namespace rrllvm return loadReactionRate(reaction); } + LLVMModelDataSymbols* modelDataSymbolsPtr = const_cast(&modelDataSymbols); + + if (modelDataSymbols.isNamedSpeciesReference(symbol)) + { + const LLVMModelDataSymbols::SpeciesReferenceInfo& info = + modelDataSymbolsPtr->getNamedSpeciesReferenceInfo(symbol); + + Value* value = mdbuilder.createStoichiometryLoad(info.row, info.column, symbol); + + if (info.type == LLVMModelDataSymbols::MultiReactantProduct) + { + std::string msg = "Mutable stochiometry for species which appear " + "multiple times in a single reaction is not currently " + "supported, species reference id: "; + msg += symbol; + throw_llvm_exception(msg); + } + + if (info.type == LLVMModelDataSymbols::Reactant) + { + // its consumed in the reaction, so has a negative in the stoich + // matrix + Value* negOne = ConstantFP::get(builder.getContext(), APFloat(-1.0)); + negOne->setName("neg_one"); + value = builder.CreateFMul(negOne, value, "neg_" + symbol); + } + + return cacheValue(symbol, args, value); + } + std::string msg = "Could not find requested symbol \'"; msg += symbol; msg += "\' in the model"; diff --git a/source/llvm/ModelResources.h b/source/llvm/ModelResources.h index 662e62bf95..98b6ed1060 100644 --- a/source/llvm/ModelResources.h +++ b/source/llvm/ModelResources.h @@ -61,6 +61,7 @@ namespace rrllvm GetEventDelayCodeGen::FunctionPtr getEventDelayPtr; EventTriggerCodeGen::FunctionPtr eventTriggerPtr; EventAssignCodeGen::FunctionPtr eventAssignPtr; + GetPiecewiseTriggerCodeGen::FunctionPtr getPiecewiseTriggerPtr; EvalVolatileStoichCodeGen::FunctionPtr evalVolatileStoichPtr; EvalConversionFactorCodeGen::FunctionPtr evalConversionFactorPtr; SetBoundarySpeciesAmountCodeGen::FunctionPtr setBoundarySpeciesAmountPtr; diff --git a/source/llvm/Random.cpp b/source/llvm/Random.cpp index 173224b81a..c5d64e0458 100644 --- a/source/llvm/Random.cpp +++ b/source/llvm/Random.cpp @@ -59,13 +59,7 @@ static Function* createGlobalMappingFunction(const char* funcName, static int64_t defaultSeed() { - int64_t seed = Config::getValue(Config::RANDOM_SEED).get(); - if (seed < 0) - { - // system time in mirsoseconds since 1970 - seed = getMicroSeconds(); - } - return seed; + return Config::getValue(Config::RANDOM_SEED).getAs(); } Random::Random(ModelGeneratorContext& ctx) @@ -585,6 +579,9 @@ void Random::setMaxTries(int maxTries) void Random::setRandomSeed(int64_t val) { + // system time in microseconds since 1970 + if (val == -1) + val = getMicroSeconds(); #ifdef RR_32BIT unsigned long maxl = std::numeric_limits::max() - 2; unsigned long seed = val % maxl; diff --git a/source/llvm/SBMLInitialValueSymbolResolver.cpp b/source/llvm/SBMLInitialValueSymbolResolver.cpp index 03103fa8b4..499eb21056 100644 --- a/source/llvm/SBMLInitialValueSymbolResolver.cpp +++ b/source/llvm/SBMLInitialValueSymbolResolver.cpp @@ -88,6 +88,37 @@ llvm::Value* SBMLInitialValueSymbolResolver::loadSymbolValue( return loadReactionRate(reaction); } + LLVMModelDataSymbols* modelDataSymbolsPtr = const_cast(&modelDataSymbols); + + if (modelDataSymbols.isNamedSpeciesReference(symbol)) + { + const LLVMModelDataSymbols::SpeciesReferenceInfo& info = + modelDataSymbolsPtr->getNamedSpeciesReferenceInfo(symbol); + + ModelDataIRBuilder mdbuilder(modelData, modelDataSymbols, builder); + Value* value = mdbuilder.createStoichiometryLoad(info.row, info.column, symbol); + + if (info.type == LLVMModelDataSymbols::MultiReactantProduct) + { + std::string msg = "Mutable stochiometry for species which appear " + "multiple times in a single reaction is not currently " + "supported, species reference id: "; + msg += symbol; + throw_llvm_exception(msg); + } + + if (info.type == LLVMModelDataSymbols::Reactant) + { + // it's consumed in the reaction, so has a negative in the stoich + // matrix + Value* negOne = ConstantFP::get(builder.getContext(), APFloat(-1.0)); + negOne->setName("neg_one"); + value = builder.CreateFMul(negOne, value, "neg_" + symbol); + } + + return cacheValue(symbol, args, value); + } + std::string msg = "Could not find requested symbol \'"; msg += symbol; msg += "\' in the model"; diff --git a/source/rrConfig.h b/source/rrConfig.h index 44c7fab930..1f6c7b5023 100644 --- a/source/rrConfig.h +++ b/source/rrConfig.h @@ -250,7 +250,7 @@ namespace rr { SIMULATEOPTIONS_MINIMUM_TIMESTEP, /** - * Specify The Maximum Time Step Size That The Internaal Integrator + * Specify The Maximum Time Step Size That The Internal Integrator * Will Use. Uses Integrator Estimated Value If <= 0. * see SimulateOptions::maximumTimeStep */ @@ -464,7 +464,6 @@ namespace rr { * for the random seed. */ RANDOM_SEED, - /** * use new numpy arrays with row/column names * experimental @@ -485,7 +484,7 @@ namespace rr { /** - * Relax SBML resrictions. + * Relax SBML restrictions. * * Allows some idiosyncrasies of e.g. JDesigner that libSBML does not * officially support. See: diff --git a/source/rrExecutableModel.h b/source/rrExecutableModel.h index b0c2564724..67d24e2324 100644 --- a/source/rrExecutableModel.h +++ b/source/rrExecutableModel.h @@ -846,12 +846,21 @@ namespace rr { */ virtual void getEventRoots(double time, const double *y, double *gdot) = 0; + virtual void getPiecewiseTriggerRoots(double time, const double* y, double* gdot) = 0; + virtual double getNextPendingEventTime(bool pop) = 0; virtual int getPendingEventSize() = 0; virtual void resetEvents() = 0; + /** + * We do root-finding for 'piecewise triggers': those times in a piecewise function that transition from one condition to the next. + */ + virtual int getNumPiecewiseTriggers() = 0; + + + /** * need a virtual destructor as object implementing this interface * can be deleted directly, i.e. diff --git a/source/rrRoadRunner.cpp b/source/rrRoadRunner.cpp index 809014d5dc..faf290dd15 100644 --- a/source/rrRoadRunner.cpp +++ b/source/rrRoadRunner.cpp @@ -128,7 +128,7 @@ namespace rr { /** * check if metabolic control analysis is valid for the model. * - * In effect, this checks that the the model is a pure + * In effect, this checks that the model is a pure * reaction-kinetics model with no rate rules, no events. * * Throws an invaid arg exception if not valid. @@ -594,8 +594,8 @@ namespace rr { RoadRunner::RoadRunner(const std::string &_compiler, const std::string &_tempDir, - const std::string &supportCodeDir) - : impl(new RoadRunnerImpl(_compiler, _tempDir, supportCodeDir)) + const std::string &supportCodeDir) + : impl(new RoadRunnerImpl(_compiler, _tempDir, supportCodeDir)) , dataVersionNumber(RR_VERSION_MAJOR * 10 + RR_VERSION_MINOR) { initLLVM(); @@ -1282,7 +1282,6 @@ namespace rr { case SelectionRecord::UNSCALED_CONTROL: dResult = getuCC(record.p1, record.p2); break; - case SelectionRecord::EIGENVALUE_REAL: { std::string species = record.p1; int index = impl->model->getFloatingSpeciesIndex(species); @@ -1347,7 +1346,7 @@ namespace rr { case SelectionRecord::INITIAL_GLOBAL_PARAMETER: impl->model->getGlobalParameterInitValues(1, &record.index, &dResult); break; - case SelectionRecord::INITIAL_COMPARTMENT: + case SelectionRecord::INITIAL_COMPARTMENT: impl->model->getCompartmentInitVolumes(1, &record.index, &dResult); break; case SelectionRecord::STOICHIOMETRY: { @@ -1362,7 +1361,6 @@ namespace rr { case SelectionRecord::TIME: dResult = getCurrentTime(); break; - default: dResult = 0.0; break; @@ -1511,7 +1509,7 @@ namespace rr { bool RoadRunner::clearModel() { // The model owns the shared library (if it exists), when the model is deleted, // its dtor unloads the shared lib. - impl->document.reset(new libsbml::SBMLDocument()); + impl->document.reset(new libsbml::SBMLDocument(3, 2)); impl->document->createModel(); if (impl->model) { impl->model = nullptr; @@ -1935,6 +1933,18 @@ namespace rr { impl->model->getStateVectorRate(impl->model->getTime(), 0); } + bool hasTime(const libsbml::ASTNode* astn) + { + if (astn->getType() == libsbml::AST_NAME_TIME) { + return true; + } + for (int c = 0; c < astn->getNumChildren(); c++) { + if (hasTime(astn->getChild(c))) { + return true; + } + } + return false; + } const ls::DoubleMatrix* RoadRunner::simulate(const SimulateOptions* opt) { get_self(); @@ -1950,6 +1960,36 @@ namespace rr { } + if (self.integrator->getName() == "gillespie") { + // Throw error if model has rate rules. + if (self.model->getNumRateRules() > 0) { + std::stringstream err; + err << "The gillespie integrator is unable to simulate a model with rate rules. Either change the rate rules into reactions, or remove them from the model. The IDs of elements with rate rules are: "; + list rrids; + self.model->getRateRuleIds(rrids); + for (auto rr = rrids.begin(); rr != rrids.end(); rr++) { + err << *rr; + if (rr != rrids.begin()) { + err << ", "; + } + } + err << "."; + throw std::domain_error(err.str()); + } + // Warn if event triggers have time in them. + const libsbml::ListOfEvents* loe = self.document->getModel()->getListOfEvents(); + for (int e = 0; e < loe->size(); e++) { + const libsbml::Event* event = loe->get(e); + const libsbml::ASTNode* trigger = event->getTrigger()->getMath(); + if (hasTime(trigger)) { + rrLog(Logger::LOG_ERROR) << "An event involving 'time' is present in this model, but time is not treated continuously in a gillespie simulation. Continuing, but the simulation will not be precise. Consider changing the trigger to involve species levels instead."; + } + if (event->isSetDelay()) { + rrLog(Logger::LOG_ERROR) << "An event with a delay is present in this model, but time is not treated continuously in a gillespie simulation. Continuing, but the delay will not be precise. Consider changing the trigger to not have a delay instead."; + } + } + } + applySimulateOptions(); const double timeStart = self.simulateOpt.start; @@ -2231,7 +2271,7 @@ namespace rr { double next = self.simulateOpt.getNext(i); double itime = self.integrator->integrate(tout, next - tout); - // the test suite is extremly sensetive to time differences, + // the test suite is extremly sensitive to time differences, // so need to use the *exact* time here. occasionally the integrator // will return a value just slightly off from the exact time // value. @@ -6900,8 +6940,11 @@ namespace rr { std::unordered_map indTolerances; // bool toleranceVector = (impl->integrator->getType("absolute_tolerance") == Setting::DOUBLEVECTOR); - Setting absTol = impl->integrator->getValue("absolute_tolerance"); - Setting::TypeId tolType = absTol.type(); + Setting absTol; + if (impl->integrator->hasValue("absolute_tolerance")) { + absTol = impl->integrator->getValue("absolute_tolerance"); + Setting::TypeId tolType = absTol.type(); + } // if absolute_tolerance is a double vector if (auto v1 = absTol.get_if>()) { @@ -7186,6 +7229,50 @@ namespace rr { } } + void RoadRunner::setSeed(long int seed, bool resetModel) { + Config::setValue(Config::RANDOM_SEED, seed); + if (resetModel) { + regenerateModel(true); + reset(SelectionRecord::TIME | + SelectionRecord::RATE | + SelectionRecord::FLOATING | + SelectionRecord::BOUNDARY | + SelectionRecord::COMPARTMENT | + SelectionRecord::GLOBAL_PARAMETER); + } + else { + impl->model->setRandomSeed(seed); + for (auto integrator : impl->integrators) { + if (integrator->getName() == "gillespie") + integrator->setValue("seed", Setting(seed)); + } + } + } + + int64_t RoadRunner::getSeed(const std::string &integratorName) { + if (integratorName.empty()) { + return impl->model->getRandomSeed(); + } + else if (integratorName == "gillespie") { + for (auto integrator : impl->integrators) { + if (integrator->getName() == integratorName) { + return integrator->getValue("seed").getAs(); + } + } + } + throw std::invalid_argument(integratorName + " is not set as the current integrator."); + } + + void RoadRunner::resetSeed() { + if (Config::getValue(Config::RANDOM_SEED).getAs() != -1) + setSeed(-1, false); + else { + for (auto integrator : impl->integrators) { + if (integrator->getName() == "gillespie") + integrator->setValue("seed", Setting(-1)); + } + } + } void writeDoubleVectorListToStream(std::ostream &out, const DoubleVectorList &results) { for (const std::vector &row: results) { diff --git a/source/rrRoadRunner.h b/source/rrRoadRunner.h index 5cfc913637..a34bc82c31 100644 --- a/source/rrRoadRunner.h +++ b/source/rrRoadRunner.h @@ -1600,6 +1600,24 @@ namespace rr { std::vector getConservedMoietyIds(); + /** + * @brief Set the value of Config::RANDOM_SEED + */ + void setSeed(long int seed, bool resetModel = true); + + /** + * @brief Returns the value of Config::RANDOM_SEED + */ + int64_t getSeed(const std::string &integratorName=""); + + + /** + * @brief Reset seed's value + * @details Set the value of Config::RANDOM_SEED to -1 + * + */ + void resetSeed(); + #ifndef SWIG // deprecated methods not SWIG'ed /** diff --git a/source/rrRoadRunnerOptions.cpp b/source/rrRoadRunnerOptions.cpp index 8e0552da06..66a34dc122 100644 --- a/source/rrRoadRunnerOptions.cpp +++ b/source/rrRoadRunnerOptions.cpp @@ -238,7 +238,7 @@ namespace rr { } } if (times.size() <= 1) { - throw std::invalid_argument("The 'times' setting must be a vector of at least two values."); + throw std::invalid_argument("The 'times' setting must be a vector of at least two values, as the first value is the time at the initial state of the model, and the second (and subsequent) times are the times the simulation progresses to."); } if (times[0] != start) { if (start == 0) //The default. @@ -255,10 +255,10 @@ namespace rr { double prev = start; for (size_t tv = 1; tv < times.size(); tv++) { double hstep = times[tv] - prev; - if (hstep <= 0) { + if (hstep < 0) { std::stringstream err; err << "The 'times' setting must be a vector of time values that start at the time value at the initial state of the model and increase along the vector. The value " - << times[tv] << " is less than or equal to the previous value of " << prev << "."; + << times[tv] << " is less than the previous value of " << prev << "."; throw std::invalid_argument(err.str()); } prev = times[tv]; diff --git a/source/rrSelectionRecord.cpp b/source/rrSelectionRecord.cpp index 1949058122..f0d9c41ded 100644 --- a/source/rrSelectionRecord.cpp +++ b/source/rrSelectionRecord.cpp @@ -325,7 +325,6 @@ rr::SelectionRecord::SelectionRecord(const std::string str) : index(-1), selectionType(UNKNOWN) { int complex; - if (is_ec(str, p1, p2)) { selectionType = ELASTICITY; diff --git a/source/rrUtils.cpp b/source/rrUtils.cpp index 95b72b13cb..00d4b6c176 100644 --- a/source/rrUtils.cpp +++ b/source/rrUtils.cpp @@ -467,7 +467,7 @@ bool hasUnimplementedTags(const std::string& descriptionFileName, const string& badtags.push_back("fbc"); badtags.push_back("FastReaction"); badtags.push_back("VolumeConcentrationRate"); - badtags.push_back("RateOf"); + //badtags.push_back("RateOf"); badtags.push_back("AssignedVariableStoichiometry"); if (integrator == "rk4" || integrator == "rk45") diff --git a/test/cxx_api_tests/ConfigTests.cpp b/test/cxx_api_tests/ConfigTests.cpp index 1a6217b44a..03631a3569 100644 --- a/test/cxx_api_tests/ConfigTests.cpp +++ b/test/cxx_api_tests/ConfigTests.cpp @@ -3,11 +3,12 @@ // #include "gtest/gtest.h" +#include "rrConfig.h" +#include "Poco/Path.h" // for Poco::Path::home, which is used in rrConfig.cpp #include #include #include -#include "rrConfig.h" -#include "Poco/Path.h" // for Poco::Path::home, which is used in rrConfig.cpp +#include using namespace rr; @@ -154,11 +155,61 @@ TEST_F(ConfigTests, SetLLJitOptLevelToAggressive){ TEST_F(ConfigTests, SetLLJitNumThreadsTo7){ // only print out the default value, which will // be different on different machines ==> bad test - std::cout << Config::getValue(Config::LLJIT_NUM_THREADS).getAs() << std::endl; Config::setValue(Config::LLJIT_NUM_THREADS, 7); int numThreads = Config::getValue(Config::LLJIT_NUM_THREADS); ASSERT_EQ(7, numThreads); } - +TEST_F(ConfigTests, SetSeedToNighMaxValues) { + uint64_t seedval = (std::numeric_limits::max)() - 5; + Config::setValue(Config::RANDOM_SEED, seedval); + rr::Setting seed_as_setting = Config::getValue(Config::RANDOM_SEED); + uint64_t seed = seed_as_setting.getAs(); + ASSERT_EQ(seedval, seed); + + ASSERT_THROW(seed_as_setting.getAs(), std::invalid_argument); + ASSERT_THROW(seed_as_setting.getAs(), std::invalid_argument); + ASSERT_THROW(seed_as_setting.getAs(), std::invalid_argument); + ASSERT_THROW(seed_as_setting.getAs(), std::invalid_argument); + string seedstr = seed_as_setting.toString(); + + seedval = (std::numeric_limits::max)() - 5; + Config::setValue(Config::RANDOM_SEED, seedval); + seed_as_setting = Config::getValue(Config::RANDOM_SEED); + seed = seed_as_setting.getAs(); + ASSERT_EQ(seedval, seed); + seed = seed_as_setting.getAs(); + ASSERT_EQ(seedval, seed); + ASSERT_THROW(seed_as_setting.getAs(), std::invalid_argument); + ASSERT_THROW(seed_as_setting.getAs(), std::invalid_argument); + ASSERT_THROW(seed_as_setting.getAs(), std::invalid_argument); + seedstr = seed_as_setting.toString(); + + seedval = (std::numeric_limits::max)() - 5; + Config::setValue(Config::RANDOM_SEED, seedval); + seed_as_setting = Config::getValue(Config::RANDOM_SEED); + seed = seed_as_setting.getAs(); + ASSERT_EQ(seedval, seed); + seed = seed_as_setting.getAs(); + ASSERT_EQ(seedval, seed); + seed = seed_as_setting.getAs(); + ASSERT_EQ(seedval, seed); + ASSERT_THROW(seed_as_setting.getAs(), std::invalid_argument); + ASSERT_THROW(seed_as_setting.getAs(), std::invalid_argument); + seedstr = seed_as_setting.toString(); + + seedval = (std::numeric_limits::max)() - 5; + Config::setValue(Config::RANDOM_SEED, seedval); + seed_as_setting = Config::getValue(Config::RANDOM_SEED); + seed = seed_as_setting.getAs(); + ASSERT_EQ(seedval, seed); + seed = seed_as_setting.getAs(); + ASSERT_EQ(seedval, seed); + seed = seed_as_setting.getAs(); + ASSERT_EQ(seedval, seed); + seed = seed_as_setting.getAs(); + ASSERT_EQ(seedval, seed); + ASSERT_THROW(seed_as_setting.getAs(), std::invalid_argument); + seedstr = seed_as_setting.toString(); +} diff --git a/test/cxx_api_tests/SettingTests.cpp b/test/cxx_api_tests/SettingTests.cpp index 9df32bd977..882bc3c5d3 100644 --- a/test/cxx_api_tests/SettingTests.cpp +++ b/test/cxx_api_tests/SettingTests.cpp @@ -201,7 +201,7 @@ TEST_F(SettingTests, ImplicitCastNegativeIntToUInt) { Setting setting(-4); ASSERT_THROW( unsigned int x = setting, - std::bad_variant_access + std::invalid_argument ); } @@ -210,7 +210,7 @@ TEST_F(SettingTests, ImplicitCastNegativeIntToULong) { Setting setting(-4); ASSERT_THROW( unsigned long x = setting, - std::bad_variant_access + std::invalid_argument ); } @@ -221,7 +221,18 @@ TEST_F(SettingTests, ImplicitCastLongToIntOutOfBounds) { Setting setting(out_of_range); ASSERT_THROW( int x = setting, - std::bad_variant_access + std::invalid_argument + ); +} + +TEST_F(SettingTests, ImplicitCastLongToNegIntOutOfBounds) { + // should raise error + std::int64_t biggest_int = std::numeric_limits::max(); + std::int64_t out_of_range = -biggest_int * 10; + Setting setting(out_of_range); + ASSERT_THROW( + int x = setting, + std::invalid_argument ); } diff --git a/test/llvm_tests/LLVMGeneratedFunctionTests.cpp b/test/llvm_tests/LLVMGeneratedFunctionTests.cpp index 38ff0b69ea..66ad2df745 100644 --- a/test/llvm_tests/LLVMGeneratedFunctionTests.cpp +++ b/test/llvm_tests/LLVMGeneratedFunctionTests.cpp @@ -51,7 +51,7 @@ void LLVMGeneratedFunctionTests::checkInitialEvaluationSpecies(int options, int // we need an instance of modelData to work with the functions. LLVMModelDataSymbols sym(m, options); std::unique_ptr r = std::make_unique(); - LLVMModelData *modelData = createModelData(sym, r.get()); + LLVMModelData *modelData = createModelData(sym, r.get(), 0); // Evaluate initial conditions evalInitConditions(modelData, options); diff --git a/test/mockups/MockExecutableModel.h b/test/mockups/MockExecutableModel.h index 71e96b38b1..5554e9bbfe 100644 --- a/test/mockups/MockExecutableModel.h +++ b/test/mockups/MockExecutableModel.h @@ -96,11 +96,13 @@ class MockExecutableModel : public rr::ExecutableModel{ MOCK_METHOD(std::string, getInfo, (), (override)); MOCK_METHOD(void, print, (std::ostream & stream), (override)); MOCK_METHOD(int, getNumEvents, (), (override)); + MOCK_METHOD(int, getNumPiecewiseTriggers, (), (override)); MOCK_METHOD(int, getEventTriggers, (size_t len, const int *indx, unsigned char *values), (override)); MOCK_METHOD(int, applyEvents, (double timeEnd, const unsigned char *previousEventStatus,const double *initialState, double *finalState), (override)); MOCK_METHOD(void, getEventRoots, (double time, const double *y, double *gdot), (override)); + MOCK_METHOD(void, getPiecewiseTriggerRoots, (double time, const double* y, double* gdot), (override)); MOCK_METHOD(double, getNextPendingEventTime, (bool pop), (override)); MOCK_METHOD(int, getPendingEventSize, (), (override)); MOCK_METHOD(void, resetEvents, (), (override)); diff --git a/test/model_analysis/model_analysis.cpp b/test/model_analysis/model_analysis.cpp index aecc99f9ea..c3310c863e 100644 --- a/test/model_analysis/model_analysis.cpp +++ b/test/model_analysis/model_analysis.cpp @@ -22,6 +22,66 @@ class ModelAnalysisTests : public RoadRunnerTest { }; +TEST_F(ModelAnalysisTests, SimulateWithSameTimes) { + RoadRunner rr((modelAnalysisModelsDir / "BIOMD0000000035_url.xml").string()); + std::vector times; + times.push_back(0); + times.push_back(1); + times.push_back(1); + times.push_back(10); + SimulateOptions opt; + opt.times = times; + const ls::DoubleMatrix* result = rr.simulate(&opt); + EXPECT_EQ(result->numRows(), 4); + EXPECT_EQ(result->Element(0, 0), 0); + EXPECT_EQ(result->Element(1, 0), 1); + EXPECT_EQ(result->Element(2, 0), 1); + EXPECT_EQ(result->Element(3, 0), 10); + for (unsigned int col = 0; col < result->numCols(); col++) { + EXPECT_EQ(result->Element(1, col), result->Element(2, col)); + } +} + +//TEST_F(ModelAnalysisTests, SimulateWithSameStartTimes) { +// RoadRunner rr((modelAnalysisModelsDir / "BIOMD0000000035_url.xml").string()); +// std::vector times; +// times.push_back(0); +// times.push_back(0); +// times.push_back(1); +// times.push_back(10); +// SimulateOptions opt; +// opt.times = times; +// const ls::DoubleMatrix* result = rr.simulate(&opt); +// EXPECT_EQ(result->numRows(), 4); +// EXPECT_EQ(result->Element(0, 0), 0); +// EXPECT_EQ(result->Element(1, 0), 0); +// EXPECT_EQ(result->Element(2, 0), 1); +// EXPECT_EQ(result->Element(3, 0), 10); +// for (unsigned int col = 0; col < result->numCols(); col++) { +// EXPECT_EQ(result->Element(0, col), result->Element(1, col)); +// } +//} + +TEST_F(ModelAnalysisTests, SimulateWithSameEndTimes) { + RoadRunner rr((modelAnalysisModelsDir / "BIOMD0000000035_url.xml").string()); + std::vector times; + times.push_back(0); + times.push_back(1); + times.push_back(10); + times.push_back(10); + SimulateOptions opt; + opt.times = times; + const ls::DoubleMatrix* result = rr.simulate(&opt); + EXPECT_EQ(result->numRows(), 4); + EXPECT_EQ(result->Element(0, 0), 0); + EXPECT_EQ(result->Element(1, 0), 1); + EXPECT_EQ(result->Element(2, 0), 10); + EXPECT_EQ(result->Element(3, 0), 10); + for (unsigned int col = 0; col < result->numCols(); col++) { + EXPECT_EQ(result->Element(2, col), result->Element(3, col)); + } +} + TEST_F(ModelAnalysisTests, issue1031) { //Config::setValue(Config::LLVM_BACKEND, Config::LLVM_BACKEND_VALUES::LLJIT); rr::RoadRunner rr((modelAnalysisModelsDir / "zero_rate_at_steady_state.xml").string()); @@ -70,9 +130,6 @@ TEST_F(ModelAnalysisTests, issue1020_full) { EXPECT_NEAR(jacobian.Element(r, c), jacobian_tv.Element(r, c), 0.0001); } } - - - } @@ -503,11 +560,11 @@ TEST_F(ModelAnalysisTests, DISABLED_SimulateCVODEFromNegativeStart_T0fireDelay) ASSERT_EQ(result->CSize(), 2); ASSERT_EQ(result->RSize(), 121); //Spot checks for when the events fire. - //EXPECT_NEAR(result->Element(0, 1), 550, 1); - //EXPECT_NEAR(result->Element(21, 1), 492.5, 1); - //EXPECT_NEAR(result->Element(22, 1), 301.7, 1); - //EXPECT_NEAR(result->Element(44, 1), 303.3, 1); - //EXPECT_NEAR(result->Element(120, 1), 400, 1); + EXPECT_NEAR(result->Element(0, 1), 550, 1); + EXPECT_NEAR(result->Element(21, 1), 492.5, 1); + EXPECT_NEAR(result->Element(22, 1), 301.7, 1); + EXPECT_NEAR(result->Element(44, 1), 303.3, 1); + EXPECT_NEAR(result->Element(120, 1), 400, 1); std::cout << *result; } @@ -516,6 +573,8 @@ TEST_F(ModelAnalysisTests, SimulateGillespieFromNegativeStart_General) { RoadRunner rr((modelAnalysisModelsDir / "negstart_event.xml").string()); rr.setIntegrator("gillespie"); rr.getIntegrator()->setValue("variable_step_size", false); + rr.getIntegrator()->setValue("seed", -1); + //std::cout << rr.getSeed("gillespie") << std::endl; SimulateOptions opt; opt.start = -2; opt.duration = 10; @@ -529,43 +588,12 @@ TEST_F(ModelAnalysisTests, SimulateGillespieFromNegativeStart_General) { } } -TEST_F(ModelAnalysisTests, SimulateGillespieFromNegativeStart_Time) { - //Event: at time > -1.1: S1 = 0; - RoadRunner rr((modelAnalysisModelsDir / "negstart_event_time.xml").string()); - rr.setIntegrator("gillespie"); - rr.getIntegrator()->setValue("variable_step_size", false); - SimulateOptions opt; - opt.start = -2; - opt.duration = 10; - opt.steps = 120; - const ls::DoubleMatrix* result = rr.simulate(&opt); - - EXPECT_LE(result->Element(11, 1), 20); - EXPECT_NEAR(result->Element(24, 1), 100, 50); -} - -TEST_F(ModelAnalysisTests, SimulateGillespieFromNegativeStart_Time_Delay) { - //Event: at 0.5 after time > -1.1: S1 = 0; - RoadRunner rr((modelAnalysisModelsDir / "negstart_event_time_delay.xml").string()); - rr.setIntegrator("gillespie"); - rr.getIntegrator()->setValue("variable_step_size", false); - SimulateOptions opt; - opt.start = -2; - opt.duration = 10; - opt.steps = 120; - const ls::DoubleMatrix* result = rr.simulate(&opt); - - EXPECT_GE(result->Element(2, 1), 500); - EXPECT_LE(result->Element(17, 1), 20); - EXPECT_NEAR(result->Element(24, 1), 66, 30); - EXPECT_NEAR(result->Element(120, 1), 946, 50); -} - TEST_F(ModelAnalysisTests, SimulateGillespieFromNegativeStart_T0fire) { //Event: at S1 > 500, t0=false: S1 = 300; RoadRunner rr((modelAnalysisModelsDir / "negstart_event_t0fire.xml").string()); rr.setIntegrator("gillespie"); rr.getIntegrator()->setValue("variable_step_size", false); + rr.getIntegrator()->setValue("seed", -1); SimulateOptions opt; opt.start = -2; opt.duration = 10; @@ -580,35 +608,6 @@ TEST_F(ModelAnalysisTests, SimulateGillespieFromNegativeStart_T0fire) { } } -TEST_F(ModelAnalysisTests, SimulateGillespieFromNegativeStart_Combo) { - //Events: E0: at 0.3 after time > -1: S1 = 0; - // E1: at time > 1: S1 = 100; - // E2: at 0.5 after S1 > 500: S1 = 300; - RoadRunner rr((modelAnalysisModelsDir / "negstart_event_combo.xml").string()); - rr.setIntegrator("gillespie"); - rr.getIntegrator()->setValue("variable_step_size", false); - SimulateOptions opt; - opt.start = -2; - opt.duration = 10; - opt.steps = 120; - const ls::DoubleMatrix* result = rr.simulate(&opt); - - EXPECT_NEAR(result->Element(2, 1), 508, 10); - EXPECT_NEAR(result->Element(8, 1), 308, 10); - EXPECT_LE(result->Element(16, 1), 15); - EXPECT_NEAR(result->Element(36, 1), 187, 20); - EXPECT_NEAR(result->Element(37, 1), 110, 10); - EXPECT_NEAR(result->Element(80, 1), 503, 30); - EXPECT_NEAR(result->Element(86, 1), 304, 10); - EXPECT_NEAR(result->Element(120, 1), 360, 40); - for (int i = 1; i <= opt.steps; i++) - { - EXPECT_LE(result->Element(i, 1), 580); - EXPECT_GE(result->Element(i, 1), 0); - } -} - - TEST_F(ModelAnalysisTests, NonZeroStarts_CVODE) { //Event: at S1 > 500, t0=false: S1 = 300; RoadRunner rr((modelAnalysisModelsDir / "BIOMD0000000035_url.xml").string()); @@ -640,8 +639,6 @@ TEST_F(ModelAnalysisTests, NonZeroStarts_CVODE) { } } - - TEST_F(ModelAnalysisTests, NonZeroStarts_RK4) { //Event: at S1 > 500, t0=false: S1 = 300; RoadRunner rr((modelAnalysisModelsDir / "threestep.xml").string()); @@ -659,7 +656,6 @@ TEST_F(ModelAnalysisTests, NonZeroStarts_RK4) { opt.start = 2; ls::DoubleMatrix s2result(*(rr.simulate(&opt))); - for (int i = 0; i <= opt.steps; i++) { EXPECT_NEAR(s0result.Element(i, 0), sneg2result.Element(i, 0) + 2, 0.01); @@ -674,7 +670,6 @@ TEST_F(ModelAnalysisTests, NonZeroStarts_RK4) { } - TEST_F(ModelAnalysisTests, NonZeroStarts_RK45) { //Event: at S1 > 500, t0=false: S1 = 300; RoadRunner rr((modelAnalysisModelsDir / "threestep.xml").string()); @@ -692,7 +687,6 @@ TEST_F(ModelAnalysisTests, NonZeroStarts_RK45) { opt.start = 2; ls::DoubleMatrix s2result(*(rr.simulate(&opt))); - for (int i = 0; i <= opt.steps; i++) { EXPECT_NEAR(s0result.Element(i, 0), sneg2result.Element(i, 0) + 2, 0.01); @@ -706,9 +700,7 @@ TEST_F(ModelAnalysisTests, NonZeroStarts_RK45) { } } - -//I don't think setting the seed is working for this one. -TEST_F(ModelAnalysisTests, DISABLED_NonZeroStarts_Gillespie) { +TEST_F(ModelAnalysisTests, NonZeroStarts_Gillespie) { RoadRunner rr((modelAnalysisModelsDir / "threestep.xml").string()); rr.setIntegrator("gillespie"); rr.getIntegrator()->setValue("variable_step_size", false); @@ -719,13 +711,14 @@ TEST_F(ModelAnalysisTests, DISABLED_NonZeroStarts_Gillespie) { opt.steps = 50; ls::DoubleMatrix s0result(*(rr.simulate(&opt))); rr.reset(int(SelectionRecord::SelectionType::ALL)); + rr.getIntegrator()->setValue("seed", 1001); opt.start = -2; ls::DoubleMatrix sneg2result(*(rr.simulate(&opt))); rr.reset(int(SelectionRecord::SelectionType::ALL)); + rr.getIntegrator()->setValue("seed", 1001); opt.start = 2; ls::DoubleMatrix s2result(*(rr.simulate(&opt))); - for (int i = 0; i <= opt.steps; i++) { EXPECT_NEAR(s0result.Element(i, 0), sneg2result.Element(i, 0) + 2, 0.01); @@ -740,7 +733,6 @@ TEST_F(ModelAnalysisTests, DISABLED_NonZeroStarts_Gillespie) { } - TEST_F(ModelAnalysisTests, GetRateOfConservedSpecies) { RoadRunner rr((modelAnalysisModelsDir / "conserved_cycle.xml").string()); rr.setConservedMoietyAnalysis(true); @@ -1699,3 +1691,58 @@ TEST_F(ModelAnalysisTests, Stoichiometry_MultiReactantProduct) { EXPECT_THROW(rr.setValue("n", 3), rrllvm::LLVMException); EXPECT_THROW(rr.setValue("stoich(S1,_J0)", 3), rrllvm::LLVMException); } + +TEST_F(ModelAnalysisTests, SetSeed_Get_Seed_Value_Check) { + RoadRunner rr((modelAnalysisModelsDir / "simple_distrib_model.xml").string()); + EXPECT_NE(rr.getSeed(), -1); + rr.setSeed(42); + EXPECT_EQ(rr.getSeed(), 42); + rr.setSeed(-1); + EXPECT_NE(rr.getSeed(), -1); + rr.setIntegrator("gillespie"); + EXPECT_NE(rr.getSeed("gillespie"), -1); + EXPECT_EQ(rr.getSeed("gillespie"), rr.getIntegrator()->getValue("seed").getAs()); + rr.setSeed(42); + EXPECT_EQ(rr.getSeed("gillespie"), 42); + EXPECT_EQ(rr.getSeed("gillespie"), rr.getIntegrator()->getValue("seed").getAs()); + rr.resetSeed(); + EXPECT_EQ(Config::getValue(Config::RANDOM_SEED).getAs(), -1); + EXPECT_NE(rr.getSeed(), -1); + EXPECT_NE(rr.getSeed("gillespie"), -1); + EXPECT_EQ(rr.getSeed("gillespie"), rr.getIntegrator()->getValue("seed").getAs()); +} + +TEST_F(ModelAnalysisTests, SetSeed_Get_Seed_Value_With_Reset_Model) { + RoadRunner rr1((modelAnalysisModelsDir / "simple_distrib_model.xml").string()); + RoadRunner rr2((modelAnalysisModelsDir / "simple_distrib_model.xml").string()); + rr1.setSeed(42); + rr2.setSeed(42); + rr1.simulate(); + rr2.simulate(); + EXPECT_NEAR(rr1.getValue("a"), rr2.getValue("a"), 0.0000001); +} + +TEST_F(ModelAnalysisTests, SetSeed_Get_Seed_Value_Uniform_Without_Reset_Model) { + RoadRunner rr1((modelAnalysisModelsDir / "uniform_distrib_model.xml").string()); + RoadRunner rr2((modelAnalysisModelsDir / "uniform_distrib_model.xml").string()); + rr1.simulate(0, 5, 50); + rr2.simulate(0, 5, 50); + rr1.setSeed(42, false); + rr2.setSeed(42, false); + rr1.simulate(5, 10, 50); + rr2.simulate(5, 10, 50); + EXPECT_NEAR(rr1.getValue("x"), rr2.getValue("x"), 0.0000001); +} + +TEST_F(ModelAnalysisTests, Set_Gillespie_Random_Seed) { + RoadRunner rr1((modelAnalysisModelsDir / "gillespie_random_seed.xml").string()); + RoadRunner rr2((modelAnalysisModelsDir / "gillespie_random_seed.xml").string()); + rr1.setIntegrator("gillespie"); + rr1.setSeed(-1, false); + rr1.simulate(0, 10, 2); + rr2.setIntegrator("gillespie"); + rr2.setSeed(-1, false); + rr2.simulate(0, 10, 2); + EXPECT_NE(rr1.getValue("S2"), rr2.getValue("S2")); +} + diff --git a/test/models/ModelAnalysis/gillespie_random_seed.xml b/test/models/ModelAnalysis/gillespie_random_seed.xml new file mode 100644 index 0000000000..d22b6b749d --- /dev/null +++ b/test/models/ModelAnalysis/gillespie_random_seed.xml @@ -0,0 +1,54 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + k1 + S1 + + + + + + + + + + + + + + + + k2 + S2 + + + + + + + diff --git a/test/models/ModelAnalysis/simple_distrib_model.xml b/test/models/ModelAnalysis/simple_distrib_model.xml new file mode 100644 index 0000000000..b70b63bf28 --- /dev/null +++ b/test/models/ModelAnalysis/simple_distrib_model.xml @@ -0,0 +1,23 @@ + + + + + + + + + + + + + + + normal + 1 + 0.5 + + + + + + diff --git a/test/models/ModelAnalysis/uniform_distrib_model.xml b/test/models/ModelAnalysis/uniform_distrib_model.xml new file mode 100644 index 0000000000..c39f260c7f --- /dev/null +++ b/test/models/ModelAnalysis/uniform_distrib_model.xml @@ -0,0 +1,46 @@ + + + + + + + + + + + + 1 + + + + + + + + + + timer + 0.99 + + + + + + + + uniform + 0 + 1 + + + + + + 0 + + + + + + + diff --git a/test/models/SBMLFeatures/need_piecewise_and_event_triggers.xml b/test/models/SBMLFeatures/need_piecewise_and_event_triggers.xml new file mode 100644 index 0000000000..47d1f583d5 --- /dev/null +++ b/test/models/SBMLFeatures/need_piecewise_and_event_triggers.xml @@ -0,0 +1,97 @@ + + + + + + + + + + + + + + + + + + + + + + 0 + + + time + 2 + + + + 3 + + + + + + + + + v + w + + + + + + + + + + + + + + + + + + k + + X1 + + u + + + + + + + + + + + + + + time + 1 + + + + time + 1.5 + + + + + + + + 3 + + + + + + + diff --git a/test/models/SBMLFeatures/need_piecewise_triggers.xml b/test/models/SBMLFeatures/need_piecewise_triggers.xml new file mode 100644 index 0000000000..6ba1ee893e --- /dev/null +++ b/test/models/SBMLFeatures/need_piecewise_triggers.xml @@ -0,0 +1,72 @@ + + + + + + + + + + + + + + + + + + + + + + + time + 1 + + + + time + 1.5 + + + + + + + + + 3 + in_interval + + + 0 + + + + + + + + + + + + + + + + + + + k + + X1 + + u + + + + + + + diff --git a/test/models/SBMLFeatures/need_piecewise_triggers_3part.xml b/test/models/SBMLFeatures/need_piecewise_triggers_3part.xml new file mode 100644 index 0000000000..42de63c60e --- /dev/null +++ b/test/models/SBMLFeatures/need_piecewise_triggers_3part.xml @@ -0,0 +1,66 @@ + + + + + + + + + + + + + + + + + + + + 0 + + + time + 1 + + + + 0 + + + time + 1.5 + + + + 3 + + + + + + + + + + + + + + + + + + + k + + X1 + + u + + + + + + + diff --git a/test/models/SBMLFeatures/need_piecewise_triggers_rev.xml b/test/models/SBMLFeatures/need_piecewise_triggers_rev.xml new file mode 100644 index 0000000000..f60b612750 --- /dev/null +++ b/test/models/SBMLFeatures/need_piecewise_triggers_rev.xml @@ -0,0 +1,66 @@ + + + + + + + + + + + + + + + + + + + + 0 + + + + + time + 1 + + + + time + 1.5 + + + + + 3 + + + + + + + + + + + + + + + + + + + k + + X1 + + u + + + + + + + diff --git a/test/models/SBMLFeatures/rateOf_1reaction.xml b/test/models/SBMLFeatures/rateOf_1reaction.xml new file mode 100644 index 0000000000..bb0123e838 --- /dev/null +++ b/test/models/SBMLFeatures/rateOf_1reaction.xml @@ -0,0 +1,37 @@ + + + + + + + + + + + + + + + + + + rateOf + s1 + + + + + + + + + + + + 1.2 + + + + + + diff --git a/test/models/SBMLFeatures/rateOf_assignmentRule1.xml b/test/models/SBMLFeatures/rateOf_assignmentRule1.xml new file mode 100644 index 0000000000..f476951098 --- /dev/null +++ b/test/models/SBMLFeatures/rateOf_assignmentRule1.xml @@ -0,0 +1,20 @@ + + + + + + + + + + + + + rateOf + s1 + + + + + + diff --git a/test/models/SBMLFeatures/rateOf_assignmentRule2.xml b/test/models/SBMLFeatures/rateOf_assignmentRule2.xml new file mode 100644 index 0000000000..fc10cf67fb --- /dev/null +++ b/test/models/SBMLFeatures/rateOf_assignmentRule2.xml @@ -0,0 +1,25 @@ + + + + + + + + + + + + 1 + + + + + + rateOf + s1 + + + + + + diff --git a/test/models/SBMLFeatures/rateOf_noARrates_invalid.xml b/test/models/SBMLFeatures/rateOf_noARrates_invalid.xml new file mode 100644 index 0000000000..e132ae4d6c --- /dev/null +++ b/test/models/SBMLFeatures/rateOf_noARrates_invalid.xml @@ -0,0 +1,27 @@ + + + + + + + + + + + + + rateOf + y + + + + + + + + 3 + + + + + diff --git a/test/models/SBMLFeatures/rateOf_noSpeciesInARComp_conc_invalid.xml b/test/models/SBMLFeatures/rateOf_noSpeciesInARComp_conc_invalid.xml new file mode 100644 index 0000000000..7b74186c56 --- /dev/null +++ b/test/models/SBMLFeatures/rateOf_noSpeciesInARComp_conc_invalid.xml @@ -0,0 +1,32 @@ + + + + + + + + + + + + + + + + + + rateOf + y + + + + + + + + 5 + + + + + diff --git a/test/models/SBMLFeatures/rateOf_recursive_assignmentRule_invalid.xml b/test/models/SBMLFeatures/rateOf_recursive_assignmentRule_invalid.xml new file mode 100644 index 0000000000..ee04ae674c --- /dev/null +++ b/test/models/SBMLFeatures/rateOf_recursive_assignmentRule_invalid.xml @@ -0,0 +1,25 @@ + + + + + + + + + + + + s2 + + + + + + rateOf + s1 + + + + + + diff --git a/test/models/SBMLFeatures/rateOf_recursive_assignmentRule_invalid2.xml b/test/models/SBMLFeatures/rateOf_recursive_assignmentRule_invalid2.xml new file mode 100644 index 0000000000..dd5a4002bb --- /dev/null +++ b/test/models/SBMLFeatures/rateOf_recursive_assignmentRule_invalid2.xml @@ -0,0 +1,28 @@ + + + + + + + + + + + + + rateOf + s2 + + + + + + + rateOf + s1 + + + + + + diff --git a/test/models/SBMLFeatures/rateOf_recursive_rateRule_invalid.xml b/test/models/SBMLFeatures/rateOf_recursive_rateRule_invalid.xml new file mode 100644 index 0000000000..6b6c51ed75 --- /dev/null +++ b/test/models/SBMLFeatures/rateOf_recursive_rateRule_invalid.xml @@ -0,0 +1,19 @@ + + + + + + + + + + + + rateOf + x + + + + + + diff --git a/test/models/SBMLFeatures/rateOf_recursive_rateRule_invalid2.xml b/test/models/SBMLFeatures/rateOf_recursive_rateRule_invalid2.xml new file mode 100644 index 0000000000..0b662f7a09 --- /dev/null +++ b/test/models/SBMLFeatures/rateOf_recursive_rateRule_invalid2.xml @@ -0,0 +1,28 @@ + + + + + + + + + + + + + rateOf + y + + + + + + + rateOf + x + + + + + + diff --git a/test/models/SBMLFeatures/rateOf_recursive_reaction_invalid.xml b/test/models/SBMLFeatures/rateOf_recursive_reaction_invalid.xml new file mode 100644 index 0000000000..8adfa522e8 --- /dev/null +++ b/test/models/SBMLFeatures/rateOf_recursive_reaction_invalid.xml @@ -0,0 +1,27 @@ + + + + + + + + + + + + + + + + + + + rateOf + x + + + + + + + diff --git a/test/models/SBMLFeatures/rateOf_recursive_reaction_invalid2.xml b/test/models/SBMLFeatures/rateOf_recursive_reaction_invalid2.xml new file mode 100644 index 0000000000..9a579209cc --- /dev/null +++ b/test/models/SBMLFeatures/rateOf_recursive_reaction_invalid2.xml @@ -0,0 +1,44 @@ + + + + + + + + + + + + + + + + + + + + + + + rateOf + y + + + + + + + + + + + + rateOf + x + + + + + + + diff --git a/test/python/test_python_api.py b/test/python/test_python_api.py index dbded7b0da..57d4a888be 100644 --- a/test/python/test_python_api.py +++ b/test/python/test_python_api.py @@ -3,6 +3,7 @@ import os import sys import unittest +import numpy thisDir = os.path.dirname(os.path.realpath(__file__)) rr_site_packages = os.path.dirname(os.path.dirname(thisDir)) @@ -1206,6 +1207,14 @@ def test_simulateWithTimes(self): result = self.rr.getSimulationData() self.assertEqual(list(result[:,0]), [0, 1, 5, 10]) + def test_simulateWithNumpyIntTimes(self): + self.rr.resetToOrigin() + numpy_intvec = numpy.array([0, 1, 5, 10]) + + self.rr.simulate(times = numpy_intvec) + result = self.rr.getSimulationData() + self.assertEqual(list(result[:,0]), [0, 1, 5, 10]) + if __name__ == "__main__": unittest.main() diff --git a/test/python/test_runStochasticTestSuite.py b/test/python/test_runStochasticTestSuite.py index 866752b435..a8db4bbcd7 100644 --- a/test/python/test_runStochasticTestSuite.py +++ b/test/python/test_runStochasticTestSuite.py @@ -71,7 +71,6 @@ def countWrongMeans(self, means, expectedmeans, expectedsds, xrange): for Z in Zvec: if Z < xrange[0] or Z > xrange[1]: nmean_wrong += 1 - self.assertLessEqual(nmean_wrong, 5) return nmean_wrong, Zvec def countWrongSDs(self, sds, expectedsds, xrange): diff --git a/test/sbml-test-suite/semantic/01822/01822-model.m b/test/sbml-test-suite/semantic/01822/01822-model.m new file mode 100644 index 0000000000..3d3ce5e373 --- /dev/null +++ b/test/sbml-test-suite/semantic/01822/01822-model.m @@ -0,0 +1,50 @@ +(* + +category: Test +synopsis: Testing the rateOf csymbol together with a non-unity changing compartment +componentTags: AssignmentRule, CSymbolRateOf, Compartment, Parameter, RateRule, Reaction, Species +testTags: Amount, InitialValueReassigned, NonConstantCompartment, NonConstantParameter, NonUnityCompartment +testType: TimeCourse +levels: 3.2 +generatedBy: Analytic +packagesPresent: comp + + In this model, a species in a non-unity compartment is affected by a reaction, is inside a changing compartment, and is set hasOnlySubstanceUnits=false, meaning that outside of reactions, it is to be treated as a concentration. This includes the 'rateOf' csymbol, meaning that 'rateOf(S1)' means the rate of change of the concentration of S1, not the rate of change of the amount of S1. Otherwise identical to test 1456, in this model, the compartment size also is changing in time. + + If you need to calculate the rate of change of the concentration, but only have the rate of change of the amount (and the rate of change of the compartment), use the following formula, where S1 is the amount, [S1] is the concentration, and C is the compartment. + + rateOf([S1]) = rateOf(S1)/C - [S1]/C * rateOf(C) + +The model contains: +* 1 species (S1) +* 1 parameter (x) +* 1 compartment (C) + +There is one reaction: + +[{width:30em,margin: 1em auto}| *Reaction* | *Rate* | +| J0: -> S1 | $2$ |] + + +There are 2 rules: + +[{width:30em,margin: 1em auto}| *Type* | *Variable* | *Formula* | +| Rate | C | $1.3$ | +| Assignment | x | $rateOf(S1)$ |] + +The initial conditions are as follows: + +[{width:35em,margin: 1em auto}| | *Value* | *Constant* | +| Initial amount of species S1 | $1$ | variable | +| Initial value of parameter x | $rateOf(S1)$ | variable | +| Initial volume of compartment 'C' | $2$ | variable |] + +The species' initial quantities are given in terms of substance units to +make it easier to use the model in a discrete stochastic simulator, but +their symbols represent their values in concentration units where they +appear in expressions. + +Note: The test data for this model was generated from an analytical +solution of the system of equations. + +*) diff --git a/test/sbml-test-suite/semantic/01822/01822-results.csv b/test/sbml-test-suite/semantic/01822/01822-results.csv new file mode 100644 index 0000000000..e4182d5a68 --- /dev/null +++ b/test/sbml-test-suite/semantic/01822/01822-results.csv @@ -0,0 +1,12 @@ +time,S1,C,x +0,1,2,0.675 +1,3,3.3,0.247933884 +2,5,4.6,0.127599244 +3,7,5.9,0.077563918 +4,9,7.2,0.052083333 +5,11,8.5,0.037370242 +6,13,9.8,0.028113286 +7,15,11.1,0.021913806 +8,17,12.4,0.017559834 +9,19,13.7,0.014385423 +10,21,15,0.012 diff --git a/test/sbml-test-suite/semantic/01822/01822-sbml-l3v2.xml b/test/sbml-test-suite/semantic/01822/01822-sbml-l3v2.xml new file mode 100644 index 0000000000..dee1e5cdec --- /dev/null +++ b/test/sbml-test-suite/semantic/01822/01822-sbml-l3v2.xml @@ -0,0 +1,42 @@ + + + + + + + + + + + + + + + + + 1.3 + + + + + + rateOf + S1 + + + + + + + + + + + + 2 + + + + + + diff --git a/test/sbml-test-suite/semantic/01822/01822-settings.txt b/test/sbml-test-suite/semantic/01822/01822-settings.txt new file mode 100644 index 0000000000..b9cd017361 --- /dev/null +++ b/test/sbml-test-suite/semantic/01822/01822-settings.txt @@ -0,0 +1,8 @@ +start: 0 +duration: 10 +steps: 10 +variables: S1, C, x +absolute: 0.0001 +relative: 0.0001 +amount: S1 +concentration: diff --git a/test/sbml-test-suite/semantic/01823/01823-model.m b/test/sbml-test-suite/semantic/01823/01823-model.m new file mode 100644 index 0000000000..62f5663e74 --- /dev/null +++ b/test/sbml-test-suite/semantic/01823/01823-model.m @@ -0,0 +1,41 @@ +(* + +category: Test +synopsis: Testing the rateOf csymbol together with a non-unity changing compartment +componentTags: AssignmentRule, CSymbolRateOf, Compartment, Parameter, RateRule, Species +testTags: Amount, InitialValueReassigned, NonConstantCompartment, NonConstantParameter, NonUnityCompartment, VolumeConcentrationRates +testType: TimeCourse +levels: 3.2 +generatedBy: Analytic +packagesPresent: comp + + In this model, a species in a non-unity changing compartment is set hasOnlySubstanceUnits=false, meaning that outside of reactions, it is to be treated as a concentration. This includes the 'rateOf' csymbol, meaning that 'rateOf(S1)' means the rate of change of the concentration of S1, not the rate of change of the amount of S1. However, it is not changing due to a reaction, but a rate rule, so its rateOf() is just the value of the rate rule. + +The model contains: +* 1 species (S1) +* 1 parameter (x) +* 1 compartment (C) + +There are 3 rules: + +[{width:30em,margin: 1em auto}| *Type* | *Variable* | *Formula* | +| Rate | S1 | $2$ | +| Rate | C | $1.3$ | +| Assignment | x | $rateOf(S1)$ |] + +The initial conditions are as follows: + +[{width:35em,margin: 1em auto}| | *Value* | *Constant* | +| Initial amount of species S1 | $1$ | variable | +| Initial value of parameter x | $rateOf(S1)$ | variable | +| Initial volume of compartment 'C' | $2$ | variable |] + +The species' initial quantities are given in terms of substance units to +make it easier to use the model in a discrete stochastic simulator, but +their symbols represent their values in concentration units where they +appear in expressions. + +Note: The test data for this model was generated from an analytical +solution of the system of equations. + +*) diff --git a/test/sbml-test-suite/semantic/01823/01823-results.csv b/test/sbml-test-suite/semantic/01823/01823-results.csv new file mode 100644 index 0000000000..f15ac19f29 --- /dev/null +++ b/test/sbml-test-suite/semantic/01823/01823-results.csv @@ -0,0 +1,12 @@ +time,S1,C,x +0,1,2,2 +1,8.25,3.3,2 +2,20.7,4.6,2 +3,38.35,5.9,2 +4,61.2,7.2,2 +5,89.25,8.5,2 +6,122.5,9.8,2 +7,160.95,11.1,2 +8,204.6,12.4,2 +9,253.45,13.7,2 +10,307.5,15,2 diff --git a/test/sbml-test-suite/semantic/01823/01823-sbml-l3v2.xml b/test/sbml-test-suite/semantic/01823/01823-sbml-l3v2.xml new file mode 100644 index 0000000000..7d9de98681 --- /dev/null +++ b/test/sbml-test-suite/semantic/01823/01823-sbml-l3v2.xml @@ -0,0 +1,35 @@ + + + + + + + + + + + + + + + + + 2 + + + + + 1.3 + + + + + + rateOf + S1 + + + + + + diff --git a/test/sbml-test-suite/semantic/01823/01823-settings.txt b/test/sbml-test-suite/semantic/01823/01823-settings.txt new file mode 100644 index 0000000000..b9cd017361 --- /dev/null +++ b/test/sbml-test-suite/semantic/01823/01823-settings.txt @@ -0,0 +1,8 @@ +start: 0 +duration: 10 +steps: 10 +variables: S1, C, x +absolute: 0.0001 +relative: 0.0001 +amount: S1 +concentration: diff --git a/test/sbml_features/CMakeLists.txt b/test/sbml_features/CMakeLists.txt index 303d02953c..32dc9997cd 100644 --- a/test/sbml_features/CMakeLists.txt +++ b/test/sbml_features/CMakeLists.txt @@ -2,10 +2,12 @@ add_test_executable( test_sbml_features test_targets ${SharedTestFiles} boolean_delay.cpp - invalid_sbml.cpp - valid_sbml.cpp events.cpp + invalid_sbml.cpp + piecewise_triggers.cpp + rateOf.cpp unsupported_packages.cpp + valid_sbml.cpp ) # make test_targets list global to all tests diff --git a/test/sbml_features/piecewise_triggers.cpp b/test/sbml_features/piecewise_triggers.cpp new file mode 100644 index 0000000000..957b0639aa --- /dev/null +++ b/test/sbml_features/piecewise_triggers.cpp @@ -0,0 +1,51 @@ +#include "gtest/gtest.h" +#include "rrRoadRunner.h" +#include "rrException.h" +#include "rrUtils.h" +#include "rrTestSuiteModelSimulation.h" +#include "sbml/SBMLTypes.h" +#include "sbml/SBMLReader.h" +#include "../test_util.h" +#include +#include "RoadRunnerTest.h" + +using namespace testing; +using namespace rr; +using namespace std; +using std::filesystem::path; + +class SBMLFeatures : public RoadRunnerTest { +public: + path SBMLFeaturesDir = rrTestModelsDir_ / "SBMLFeatures"; + SBMLFeatures() = default; +}; + +TEST_F(SBMLFeatures, PIECEWISE_TRIGGERS_1) +{ + RoadRunner rri((SBMLFeaturesDir / "need_piecewise_triggers.xml").string()); + const ls::DoubleMatrix* results = rri.simulate(); + EXPECT_NEAR(rri.getValue(rri.createSelection("X1")), 1.0310384266045054, 0.0001); +} + +TEST_F(SBMLFeatures, PIECEWISE_TRIGGERS_2) +{ + RoadRunner rri((SBMLFeaturesDir / "need_piecewise_triggers_rev.xml").string()); + const ls::DoubleMatrix* results = rri.simulate(); + EXPECT_NEAR(rri.getValue(rri.createSelection("X1")), 1.0310384266045054, 0.0001); +} + +TEST_F(SBMLFeatures, PIECEWISE_TRIGGERS_3) +{ + RoadRunner rri((SBMLFeaturesDir / "need_piecewise_triggers_3part.xml").string()); + const ls::DoubleMatrix* results = rri.simulate(); + EXPECT_NEAR(rri.getValue(rri.createSelection("X1")), 1.0310384266045054, 0.0001); +} + +TEST_F(SBMLFeatures, PIECEWISE_AND_EVENT_TRIGGERS) +{ + //This worked in rr 2.5, but failed with the new version because piecewise triggers overrode the events triggers. + RoadRunner rri((SBMLFeaturesDir / "need_piecewise_and_event_triggers.xml").string()); + const ls::DoubleMatrix* results = rri.simulate(); + EXPECT_NEAR(rri.getValue(rri.createSelection("X1")), 17.66586981819867, 0.0001); +} + diff --git a/test/sbml_features/rateOf.cpp b/test/sbml_features/rateOf.cpp new file mode 100644 index 0000000000..b7b8a6c0ba --- /dev/null +++ b/test/sbml_features/rateOf.cpp @@ -0,0 +1,182 @@ +#include "gtest/gtest.h" +#include "rrRoadRunner.h" +#include "rrException.h" +#include "rrUtils.h" +#include "rrTestSuiteModelSimulation.h" +#include "sbml/SBMLTypes.h" +#include "sbml/SBMLReader.h" +#include "../test_util.h" +#include +#include "RoadRunnerTest.h" + +using namespace testing; +using namespace rr; +using namespace std; +using std::filesystem::path; + +class SBMLFeatures : public RoadRunnerTest { +public: + path SBMLFeaturesDir = rrTestModelsDir_ / "SBMLFeatures"; + SBMLFeatures() = default; +}; + +TEST_F(SBMLFeatures, RATEOF_RNX1) +{ + //Logger::enableConsoleLogging(); + //Logger::setLevel(Logger::LOG_DEBUG); + + try + { + RoadRunner rri((SBMLFeaturesDir / "rateOf_1reaction.xml").string()); + rri.getSimulateOptions().duration = 2; + rri.simulate(); + EXPECT_EQ(rri.getValue(rri.createSelection("s2")), 1.2); + } + catch (std::exception& ex) + { + std::cout << "Exception: " << ex.what() << std::endl; + EXPECT_TRUE(false); + } +} + +TEST_F(SBMLFeatures, RATEOF_AR1) +{ + try + { + RoadRunner rri((SBMLFeaturesDir / "rateOf_assignmentRule1.xml").string()); + rri.getSimulateOptions().duration = 2; + rri.simulate(); + EXPECT_EQ(rri.getValue(rri.createSelection("s2")), 0.0); + } + catch (std::exception& ex) + { + std::cout << "Exception: " << ex.what() << std::endl; + EXPECT_TRUE(false); + } +} + +TEST_F(SBMLFeatures, RATEOF_AR2) +{ + try + { + RoadRunner rri((SBMLFeaturesDir / "rateOf_assignmentRule2.xml").string()); + EXPECT_EQ(rri.getValue(rri.createSelection("s2")), 1.0); + EXPECT_EQ(rri.getValue(rri.createSelection("s1'")), 1.0); + rri.getSimulateOptions().duration = 1; + rri.simulate(); + EXPECT_EQ(rri.getValue(rri.createSelection("s2")), 1.0); + EXPECT_NEAR(rri.getValue(rri.createSelection("s1")), 6.0, 0.0001); + } + catch (std::exception& ex) + { + std::cout << "Exception: " << ex.what() << std::endl; + EXPECT_TRUE(false); + } +} + +TEST_F(SBMLFeatures, RATEOF_AR_recursive_err) +{ + try + { + RoadRunner rri((SBMLFeaturesDir / "rateOf_recursive_assignmentRule_invalid.xml").string()); + EXPECT_TRUE(false); + } + catch (std::exception& ex) + { + EXPECT_TRUE(true); + } +} + +TEST_F(SBMLFeatures, RATEOF_AR_recursive_err2) +{ + try + { + RoadRunner rri((SBMLFeaturesDir / "rateOf_recursive_assignmentRule_invalid2.xml").string()); + EXPECT_TRUE(false); + } + catch (std::exception& ex) + { + EXPECT_TRUE(true); + } +} + +TEST_F(SBMLFeatures, RATEOF_RateRule_recursive_err) +{ + try + { + RoadRunner rri((SBMLFeaturesDir / "rateOf_recursive_rateRule_invalid.xml").string()); + EXPECT_TRUE(false); + } + catch (std::exception& ex) + { + EXPECT_TRUE(true); + } +} + +TEST_F(SBMLFeatures, RATEOF_RateRule_recursive_err2) +{ + try + { + RoadRunner rri((SBMLFeaturesDir / "rateOf_recursive_rateRule_invalid2.xml").string()); + EXPECT_TRUE(false); + } + catch (std::exception& ex) + { + EXPECT_TRUE(true); + } +} + +TEST_F(SBMLFeatures, RATEOF_Reaction_recursive_err) +{ + try + { + RoadRunner rri((SBMLFeaturesDir / "rateOf_recursive_reaction_invalid.xml").string()); + EXPECT_TRUE(false); + } + catch (std::exception& ex) + { + EXPECT_TRUE(true); + } +} + +TEST_F(SBMLFeatures, RATEOF_Reaction_recursive_err2) +{ + try + { + RoadRunner rri((SBMLFeaturesDir / "rateOf_recursive_reaction_invalid2.xml").string()); + EXPECT_TRUE(false); + } + catch (std::exception& ex) + { + EXPECT_TRUE(true); + } +} + +TEST_F(SBMLFeatures, RATEOF_NoRateOfAssignmentRule_err) +{ + try + { + RoadRunner rri((SBMLFeaturesDir / "rateOf_noARrates_invalid.xml").string()); + EXPECT_TRUE(false); + } + catch (std::exception& ex) + { + std::cout << "Exception: " << ex.what() << std::endl; + EXPECT_TRUE(true); + } +} + +TEST_F(SBMLFeatures, RATEOF_NoRateOfConcentrationWhenCompartmentHasAR_err) +{ + try + { + RoadRunner rri((SBMLFeaturesDir / "rateOf_noSpeciesInARComp_conc_invalid.xml").string()); + EXPECT_TRUE(false); + } + catch (std::exception& ex) + { + std::cout << "Exception: " << ex.what() << std::endl; + EXPECT_TRUE(true); + } +} + diff --git a/test/semantic_STS/sbml_test_suite.cpp b/test/semantic_STS/sbml_test_suite.cpp index c0ebbfcdf2..e73898908a 100644 --- a/test/semantic_STS/sbml_test_suite.cpp +++ b/test/semantic_STS/sbml_test_suite.cpp @@ -252,7 +252,7 @@ class SbmlTestSuite : public RoadRunnerTest TEST_F(SbmlTestSuite, test_single) { // Use when need to run one test. - EXPECT_TRUE(RunTest(1511)); + EXPECT_TRUE(RunTest(1823)); } TEST_F(SbmlTestSuite, t1) { @@ -7556,15 +7556,15 @@ TEST_F(SbmlTestSuite, t1821) { EXPECT_TRUE(RunTest(1821)); } +TEST_F(SbmlTestSuite, t1822) +{ + EXPECT_TRUE(RunTest(1822)); +} +TEST_F(SbmlTestSuite, t1823) +{ + EXPECT_TRUE(RunTest(1823)); +} //The following tests do not yet exist. -//TEST_F(SbmlTestSuite, t1822) -//{ -// EXPECT_TRUE(RunTest(1822)); -//} -//TEST_F(SbmlTestSuite, t1823) -//{ -// EXPECT_TRUE(RunTest(1823)); -//} //TEST_F(SbmlTestSuite, t1824) //{ // EXPECT_TRUE(RunTest(1824)); diff --git a/wrappers/C/rrc_api.cpp b/wrappers/C/rrc_api.cpp index 831df0c844..91b24e9df6 100644 --- a/wrappers/C/rrc_api.cpp +++ b/wrappers/C/rrc_api.cpp @@ -3410,37 +3410,6 @@ vector rr_getRatesOfChange(RoadRunner* rr) return result; } -C_DECL_SPEC bool rrcCallConv getSeed(RRHandle h, long* result) { - start_try - RoadRunner *r = (RoadRunner*)h; - //Integrator *intg = r->getIntegrator(Integrator::GILLESPIE); - //*result = intg->getItem("seed").get(); - *result = r->getIntegrator()->getValue("seed").get(); - return true; - catch_bool_macro; -} - -C_DECL_SPEC bool rrcCallConv setSeed(RRHandle h, long result) { - start_try - RoadRunner *r = (RoadRunner*)h; - //Integrator *intg = r->getIntegrator(Integrator::GILLESPIE); - //intg->setItem("seed", result); - Integrator *intg = r->getIntegrator(); - if (intg->getName() == "gillespie") - { - intg->setValue("seed", Setting(result)); - } - else - { - Integrator *intg = dynamic_cast( - IntegratorFactory::getInstance().New("gillespie", r->getModel()) - ); - intg->setValue("seed", Setting(result)); - } - return true; - catch_bool_macro -} - C_DECL_SPEC RRCDataPtr rrcCallConv gillespie(RRHandle handle) { start_try RoadRunner *r = (RoadRunner*)handle; @@ -3646,6 +3615,7 @@ bool rrcCallConv resetAll(RRHandle handle) SelectionRecord::BOUNDARY | SelectionRecord::COMPARTMENT | SelectionRecord::GLOBAL_PARAMETER); + rri->resetSeed(); return true; catch_bool_macro } diff --git a/wrappers/C/rrc_api.h b/wrappers/C/rrc_api.h index e64114e8bf..f8b1dda41e 100644 --- a/wrappers/C/rrc_api.h +++ b/wrappers/C/rrc_api.h @@ -2923,26 +2923,6 @@ C_DECL_SPEC bool rrcCallConv getuEE(RRHandle handle, const char* name, const cha // Stochastic methods // -------------------------------------------------------------------------------- -/*! - \brief Determine the current seed used by the random generator - - \param[in] handle Handle to a RoadRunner instance - \param[out] seed This is the value of the current seed, returned to the caller - \return Returns true if successful - \ingroup stochastic -*/ -C_DECL_SPEC bool rrcCallConv getSeed(RRHandle handle, long* seed); - -/*! - \brief Set the current seed used by the random generator - - \param[in] handle Handle to a RoadRunner instance - \param[out] seed This is the value the caller will set the seed to - \return Returns true if successful - \ingroup stochastic -*/ -C_DECL_SPEC bool rrcCallConv setSeed(RRHandle handle, long seed); - /*! \brief Carry out a time-course simulation using the Gillespie algorithm with variable step size. setTimeStart, setTimeEnd, etc are used to set the simulation @@ -3250,7 +3230,7 @@ C_DECL_SPEC RRStringArrayPtr rrcCallConv getListOfConfigKeys(); /*! \brief Set the selection list for output from simulate(void) or simulateEx(void) - Use setTimeCourseSelectionListEx(handle, length, list) to set the the simulate selection list. + Use setTimeCourseSelectionListEx(handle, length, list) to set the simulate selection list. Compared to setTimeCourseSelectionList, setTimeCourseSelectionListEx, expects a list of char* strings otherwise it has identical functionality. \param[in] handle Handle to a RoadRunner instance diff --git a/wrappers/C/rrc_utilities.h b/wrappers/C/rrc_utilities.h index 6ecc70f74e..249cfdb92c 100644 --- a/wrappers/C/rrc_utilities.h +++ b/wrappers/C/rrc_utilities.h @@ -76,7 +76,7 @@ extern const char* ALLOCATE_API_ERROR_MSG; extern const char* INVALID_HANDLE_ERROR_MSG; /*! - \brief Retrieves the the content of a file. + \brief Retrieves the content of a file. \param[in] fName Pointer to a string holding the name and optionla path to a file \return Content of file as a string, returns null if it fails \ingroup utilities diff --git a/wrappers/Python/roadrunner/roadrunner.i b/wrappers/Python/roadrunner/roadrunner.i index 390eb81937..1f80f97a0e 100644 --- a/wrappers/Python/roadrunner/roadrunner.i +++ b/wrappers/Python/roadrunner/roadrunner.i @@ -1203,7 +1203,7 @@ namespace std { class ostream{}; } if (PyFloat_Check(pval)) { val = PyFloat_AsDouble(pval); } - else if (PyInt_Check(pval)) { + else if (PyArray_IsIntegerScalar(pval)) { val = PyInt_AsLong(pval); } else { @@ -1281,13 +1281,48 @@ namespace std { class ostream{}; } else: return _roadrunner.RoadRunner__getValue(self, *args) - def setValues(self, keys, values): - for key, val in zip(keys, values): - _roadrunner.RoadRunner_setValue(self, key, val) + def setValues(self, keysOrDict, values=None): + """ + Sets a number of values in the roadrunner object all at once. Use either with the first argument defined as a dictionary, or with both arguments defined, with the first as the keys and the second as the values. + :param keysOrDict (list or dict): Either a list of id strings to set, or a dictionary with string keys and numerical values. + :param values (list): The list of values to use. Must be identical in length to 'keysOrDict', and keysOrDict must not be a dictionary. + """ + if isinstance(keysOrDict, dict): + if values is not None: + raise ValueError("Because keysOrDict is a dictionary, 'values' must be None.") + for key in keysOrDict: + _roadrunner.RoadRunner_setValue(self, key, keysOrDict[key]) + else: + for key, val in zip(keysOrDict, values): + _roadrunner.RoadRunner_setValue(self, key, val) def getModel(self): return self._getModel() + def setSeed(self, seed, resetModel=True): + """ + Sets the seed for all random number generators: any use of 'distrib' in the model, any + simultaneously-firing events, and any use of the 'gillespie' integrator. Also sets + the global seed in the configuration object, so subsequently-created roadrunner objects + will be created with the same seed. Setting the seed to -1 (or any negative integer) + tells the random number generators to use a seed based on the system clock (in + microcseconds) instead. + warnings.warn("the integrator option is now ignored for this function. So this function now sets the seed used for\ + the existing model and for the global configuration option") + :param seed: The seed to use for the random number generator. + :param resetModel: If True, the model will be reset after setting the seed. + """ + _roadrunner.RoadRunner_setSeed(self, seed, resetModel) + + def getSeed(self, integratorName=""): + """ + Obtain the current seed for this roadrunner object. + warnings.warn("the integrator option is now ignored for this function. So this function now returns the seed used for\ + the existing model and for the global configuration option") + :param seed: The integrator name to get the seed for. If None, the Config seed is returned. + """ + return _roadrunner.RoadRunner_getSeed(self, integratorName) + def _setConservedMoietyAnalysisProxy(self, value): self._setConservedMoietyAnalysis(value) self._makeProperties() @@ -1648,7 +1683,7 @@ namespace std { class ostream{}; } DEPRECATED: use solver wrappers (this setting only available for some solvers). seed - DEPRECATED: use solver wrappers (this setting only available for some solvers). + DEPRECATED: use 'setSeed'. :returns: a numpy array with each selected output time series being a