Let froot depend on w

WIP

This should hopefully allow for more efficient computation of `root`-derivatives
by avoiding flattening `w` into `root` which, so far, made computing `drootdt_total`
prohibitively expensive in case of large `w` dependencies in `root`.

Author: dweindl

cmake/AmiciFindBLAS.cmake

@@ -15,6 +15,11 @@ if(DEFINED ENV{AMICI_BLAS_USE_SCIPY_OPENBLAS})
       "Using AMICI_BLAS_USE_SCIPY_OPENBLAS=${AMICI_BLAS_USE_SCIPY_OPENBLAS} from environment variable."
   )
   set(AMICI_BLAS_USE_SCIPY_OPENBLAS $ENV{AMICI_BLAS_USE_SCIPY_OPENBLAS})
+elseif(NOT DEFINED AMICI_BLAS_USE_SCIPY_OPENBLAS
+       AND NOT AMICI_PYTHON_BUILD_EXT_ONLY)
+  # If were are not building the Python extension, it's unlikely that we want to
+  # use scipy-openblas
+  set(AMICI_BLAS_USE_SCIPY_OPENBLAS FALSE)
 endif()
 
 if((${BLAS} STREQUAL "MKL" OR DEFINED ENV{MKLROOT})

include/amici/abstract_model.h

@@ -299,6 +299,7 @@ class AbstractModel {
      * @param p parameter vector
      * @param k constant vector
      * @param h Heaviside vector
+     * @param w vector with helper variables
      * @param dx time derivative of state (DAE only)
      * @param tcl total abundances for conservation laws
      * @param sx current state sensitivity
@@ -307,8 +308,9 @@ class AbstractModel {
      */
     virtual void fstau(
         realtype* stau, realtype t, realtype const* x, realtype const* p,
-        realtype const* k, realtype const* h, realtype const* dx,
-        realtype const* tcl, realtype const* sx, int ip, int ie
+        realtype const* k, realtype const* h, realtype const* w,
+        realtype const* dx, realtype const* tcl, realtype const* sx, int ip,
+        int ie
     );
 
     /**
@@ -542,6 +544,7 @@ class AbstractModel {
      * @param p parameter vector
      * @param k constant vector
      * @param h Heaviside vector
+     * @param w vector with helper variables
      * @param dx time derivative of state (DAE only)
      * @param ie event index
      * @param xdot new model right hand side
@@ -552,9 +555,10 @@ class AbstractModel {
      */
     virtual void fdeltaxB(
         realtype* deltaxB, realtype t, realtype const* x, realtype const* p,
-        realtype const* k, realtype const* h, realtype const* dx, int ie,
-        realtype const* xdot, realtype const* xdot_old, realtype const* x_old,
-        realtype const* xB, realtype const* tcl
+        realtype const* k, realtype const* h, realtype const* w,
+        realtype const* dx, int ie, realtype const* xdot,
+        realtype const* xdot_old, realtype const* x_old, realtype const* xB,
+        realtype const* tcl
     );
 
     /**
@@ -565,6 +569,7 @@ class AbstractModel {
      * @param p parameter vector
      * @param k constant vector
      * @param h Heaviside vector
+     * @param w vector with helper variables
      * @param dx time derivative of state (DAE only)
      * @param ip sensitivity index
      * @param ie event index
@@ -575,9 +580,9 @@ class AbstractModel {
      */
     virtual void fdeltaqB(
         realtype* deltaqB, realtype t, realtype const* x, realtype const* p,
-        realtype const* k, realtype const* h, realtype const* dx, int ip,
-        int ie, realtype const* xdot, realtype const* xdot_old,
-        realtype const* x_old, realtype const* xB
+        realtype const* k, realtype const* h, realtype const* w,
+        realtype const* dx, int ip, int ie, realtype const* xdot,
+        realtype const* xdot_old, realtype const* x_old, realtype const* xB
     );
 
     /**

include/amici/model_dae.h

@@ -359,11 +359,13 @@ class Model_DAE : public Model {
      * @param p parameter vector
      * @param k constants vector
      * @param h Heaviside vector
+     * @param w vector with helper variables
      * @param dx Vector with the derivative states
      **/
     virtual void froot(
         realtype* root, realtype t, realtype const* x, double const* p,
-        double const* k, realtype const* h, realtype const* dx
+        double const* k, realtype const* h, realtype const* w,
+        realtype const* dx
     );
 
     /**

include/amici/model_ode.h

@@ -316,11 +316,13 @@ class Model_ODE : public Model {
      * @param p parameter vector
      * @param k constants vector
      * @param h Heaviside vector
+     * @param w vector with helper variables
      * @param tcl total abundances for conservation laws
      **/
     virtual void froot(
         realtype* root, realtype t, realtype const* x, realtype const* p,
-        realtype const* k, realtype const* h, realtype const* tcl
+        realtype const* k, realtype const* h, realtype const* w,
+        realtype const* tcl
     );
 
     /**

models/model_calvetti_py/CMakeLists.txt

@@ -4,6 +4,12 @@ cmake_policy(VERSION 3.22...3.31)
 
 project(model_calvetti_py)
 
+message(STATUS "CMAKE_CURRENT_SOURCE_DIR: ${CMAKE_CURRENT_SOURCE_DIR}")
+message(STATUS "CMAKE_BINARY_DIR: ${CMAKE_BINARY_DIR}")
+message(STATUS "CMAKE_INSTALL_PREFIX: ${CMAKE_INSTALL_PREFIX}")
+message(STATUS "CMAKE_VERSION: ${CMAKE_VERSION}")
+message(STATUS "CMAKE_COMMAND: ${CMAKE_COMMAND}")
+
 set(CMAKE_CXX_STANDARD 20)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 set(CMAKE_POSITION_INDEPENDENT_CODE ON)

models/model_calvetti_py/model_calvetti_py.h

@@ -38,7 +38,7 @@ extern void dJydy_rowvals_model_calvetti_py(SUNMatrixWrapper &rowvals, int index
 
 
 
-extern void root_model_calvetti_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl);
+extern void root_model_calvetti_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl);
 
 
 
@@ -200,10 +200,10 @@ class Model_model_calvetti_py : public amici::Model_DAE {
     void fdeltasx(realtype *deltasx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *sx, const realtype *stau, const realtype *tcl, const realtype *x_old) override {}
 
 
-    void fdeltaxB(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl) override {}
+    void fdeltaxB(realtype *deltaxB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB, const realtype *tcl) override {}
 
 
-    void fdeltaqB(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB) override {}
+    void fdeltaqB(realtype *deltaqB, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const int ip, const int ie, const realtype *xdot, const realtype *xdot_old, const realtype *x_old, const realtype *xB) override {}
 
 
     void fdrzdp(realtype *drzdp, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const int ip) override {}
@@ -323,8 +323,8 @@ class Model_model_calvetti_py : public amici::Model_DAE {
     void fdzdx(realtype *dzdx, const int ie, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h) override {}
 
 
-    void froot(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl) override {
-        root_model_calvetti_py(root, t, x, p, k, h, tcl);
+    void froot(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl) override {
+        root_model_calvetti_py(root, t, x, p, k, h, w, tcl);
     }
 
 
@@ -339,7 +339,7 @@ class Model_model_calvetti_py : public amici::Model_DAE {
     void fsigmaz(realtype *sigmaz, const realtype t, const realtype *p, const realtype *k) override {}
 
 
-    void fstau(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie) override {}
+    void fstau(realtype *stau, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *dx, const realtype *tcl, const realtype *sx, const int ip, const int ie) override {}
 
     void fsx0(realtype *sx0, const realtype t, const realtype *x, const realtype *p, const realtype *k, const int ip) override {}
 
@@ -557,7 +557,7 @@ class Model_model_calvetti_py : public amici::Model_DAE {
      * @return AMICI git commit hash
      */
     std::string get_amici_commit() const override {
-        return "f005fac9e2de7c3c90be2ac55d4ad165471ed1e7";
+        return "unknown";
     }
 
     bool has_quadratic_llh() const override {

models/model_calvetti_py/root.cpp

@@ -5,11 +5,12 @@
 #include "x.h"
 #include "k.h"
 #include "h.h"
+#include "w.h"
 
 namespace amici {
 namespace model_model_calvetti_py {
 
-void root_model_calvetti_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *tcl){
+void root_model_calvetti_py(realtype *root, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl){
     root[0] = t - 10;
     root[1] = 10 - t;
     root[2] = 12 - t;

models/model_calvetti_py/setup.py

@@ -1,5 +1,6 @@
 """AMICI model package setup"""
 
+import importlib.metadata
 import os
 import sys
 from pathlib import Path
@@ -8,7 +9,6 @@
 from amici.custom_commands import AmiciBuildCMakeExtension
 from cmake_build_extension import CMakeExtension
 from setuptools import find_namespace_packages, setup
-import importlib.metadata
 
 
 def get_extension() -> CMakeExtension:

models/model_calvetti_py/swig/model_calvetti_py.i

@@ -8,10 +8,11 @@ import sysconfig
 from pathlib import Path
 
 ext_suffix = sysconfig.get_config_var('EXT_SUFFIX')
+extension_path = Path(__file__).parent / f'_model_calvetti_py{ext_suffix}'
 _model_calvetti_py = amici._module_from_path(
     'model_calvetti_py._model_calvetti_py' if __package__ or '.' in __name__
     else '_model_calvetti_py',
-    Path(__file__).parent / f'_model_calvetti_py{ext_suffix}',
+    extension_path,
 )
 
 def _get_import_time():
@@ -36,6 +37,28 @@ if t_imported < t_modified:
 
 %module(package="model_calvetti_py",moduleimport=MODULEIMPORT) model_calvetti_py
 
+// store swig version
+%constant int SWIG_VERSION_MAJOR = (SWIG_VERSION >> 16);
+%constant int SWIG_VERSION_MINOR = ((SWIG_VERSION >> 8) & 0xff);
+%constant int SWIG_VERSION_PATCH = (SWIG_VERSION & 0xff);
+
+%pythoncode %{
+# SWIG version used to build the model extension as `(major, minor, patch)`
+_SWIG_VERSION = (SWIG_VERSION_MAJOR, SWIG_VERSION_MINOR, SWIG_VERSION_PATCH)
+
+if (amici_swig := amici.amici._SWIG_VERSION) != (model_swig := _SWIG_VERSION):
+    import warnings
+    warnings.warn(
+        f"SWIG version mismatch between amici ({amici_swig}) and model "
+        f"({model_swig}). This may lead to unexpected behavior. "
+        "In that case, please recompile the model with swig=="
+        f"{amici_swig[0]}.{amici_swig[1]}.{amici_swig[2]} or rebuild amici "
+        f"with swig=={model_swig[0]}.{model_swig[1]}.{model_swig[2]}.",
+        RuntimeWarning,
+        stacklevel=2,
+    )
+%}
+
 %pythoncode %{
 # the model-package __init__.py module (will be set during import)
 _model_module = None
@@ -56,7 +79,7 @@ using namespace amici;
 // store the time a module was imported
 %{
 #include <chrono>
-static std::chrono::time_point<std::chrono::system_clock> _module_import_time;
+static std::chrono::time_point<std::chrono::system_clock> _module_import_time = std::chrono::system_clock::now();
 
 static double _get_import_time() {
     auto epoch = _module_import_time.time_since_epoch();
@@ -67,7 +90,9 @@ static double _get_import_time() {
 static double _get_import_time();
 
 %init %{
-    _module_import_time = std::chrono::system_clock::now();
+    // NOTE: from SWIG 4.4.0 onwards, %init code is executed every time the
+    // module is executed - not only on first import
+    // This code ends up in `SWIG_mod_exec`.
 %}
 
 

models/model_dirac_py/CMakeLists.txt

@@ -4,6 +4,12 @@ cmake_policy(VERSION 3.22...3.31)
 
 project(model_dirac_py)
 
+message(STATUS "CMAKE_CURRENT_SOURCE_DIR: ${CMAKE_CURRENT_SOURCE_DIR}")
+message(STATUS "CMAKE_BINARY_DIR: ${CMAKE_BINARY_DIR}")
+message(STATUS "CMAKE_INSTALL_PREFIX: ${CMAKE_INSTALL_PREFIX}")
+message(STATUS "CMAKE_VERSION: ${CMAKE_VERSION}")
+message(STATUS "CMAKE_COMMAND: ${CMAKE_COMMAND}")
+
 set(CMAKE_CXX_STANDARD 20)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 set(CMAKE_POSITION_INDEPENDENT_CODE ON)