diff --git a/.github/workflows/sbml_test_suite.yml b/.github/workflows/sbml_test_suite.yml index c3d0a7a4..f34c5b19 100644 --- a/.github/workflows/sbml_test_suite.yml +++ b/.github/workflows/sbml_test_suite.yml @@ -20,18 +20,18 @@ jobs: strategy: fail-fast: false - # matrix is split by number of passing test cases (~230 per batch) + # matrix is split by number of passing test cases (269 per batch) matrix: include: - first_case: "1" - last_case: "245" - - first_case: "246" - last_case: "578" - - first_case: "579" - last_case: "890" - - first_case: "891" - last_case: "1267" - - first_case: "1268" + last_case: "276" + - first_case: "277" + last_case: "678" + - first_case: "679" + last_case: "1030" + - first_case: "1031" + last_case: "1364" + - first_case: "1365" last_case: "1821" steps: diff --git a/chaste_sbml/SbmlRefModels/CMakeLists.txt b/chaste_sbml/SbmlRefModels/CMakeLists.txt index 9ae069b5..88a617a2 100644 --- a/chaste_sbml/SbmlRefModels/CMakeLists.txt +++ b/chaste_sbml/SbmlRefModels/CMakeLists.txt @@ -1,2 +1,10 @@ find_package(Chaste COMPONENTS cell_based) chaste_do_project(SbmlRefModels) + +# Precompile the heavy header every generated case model includes, to speed up rebuilds. The +# serialization export wrappers are deliberately excluded (they are order-sensitive and belong in +# each .cpp, not a shared PCH). +if (TARGET chaste_project_SbmlRefModels) + target_precompile_headers(chaste_project_SbmlRefModels PRIVATE + "${CMAKE_CURRENT_SOURCE_DIR}/src/AbstractSbmlOdeSystem.hpp") +endif() diff --git a/chaste_sbml/SbmlRefModels/generate_cases.py b/chaste_sbml/SbmlRefModels/generate_cases.py index 044c052b..e20e57e0 100644 --- a/chaste_sbml/SbmlRefModels/generate_cases.py +++ b/chaste_sbml/SbmlRefModels/generate_cases.py @@ -9,7 +9,7 @@ from typing import TYPE_CHECKING from chaste_sbml import ChasteSbmlModel -from chaste_sbml._config import ROOT_DIR, ModelType +from chaste_sbml._config import ModelType from chaste_sbml._names import generate_header_guard if TYPE_CHECKING: @@ -20,6 +20,24 @@ sbml_versions = ["l2v5", "l3v2", "l3v1"] +def _cpp_double_literal(value: str) -> str: + """Render a reference-CSV cell as a C++ double literal. + + Reference results can be non-finite -- a division by zero gives INF and 0/0 gives NaN -- and + the SBML Test Suite writes these as ``INF``/``-INF``/``NaN``, which are not valid C++ tokens. + Map them to the ```` macros so the expected-results initialiser list compiles; finite + numbers are already valid literals and pass through unchanged. + """ + token = value.strip() + magnitude = token.lstrip("+-").lower() + sign = "-" if token.startswith("-") else "" + if magnitude == "inf": + return f"{sign}INFINITY" + if magnitude == "nan": + return "NAN" + return token + + class TestType(Enum): """Enumeration of test types for code generation.""" @@ -50,7 +68,9 @@ def __init__(self, sbml_file: str, sbml_version: str, test_params: dict[str, "An self._test_hpp_filename = f"Test{self._model_name}.hpp" test_result_columns = ", ".join(f'"{col}"' for col in test_params["results"][0]) - test_result_data = ["{ " + ", ".join(row) + " }" for row in test_params["results"][1:]] + test_result_data = [ + "{ " + ", ".join(_cpp_double_literal(cell) for cell in row) + " }" for row in test_params["results"][1:] + ] test_result_data = ",\n".join(test_result_data) test_amounts = test_params["settings"]["amount"].split(",") diff --git a/chaste_sbml/SbmlRefModels/infra/configure.sh b/chaste_sbml/SbmlRefModels/infra/configure.sh index 5f0dd878..485447ca 100755 --- a/chaste_sbml/SbmlRefModels/infra/configure.sh +++ b/chaste_sbml/SbmlRefModels/infra/configure.sh @@ -19,7 +19,7 @@ fi if [[ ! -d "${chaste_dir}" ]]; then echo "Error: Chaste source directory not found at '${chaste_dir}'." >&2 - echo "Set CHASTE_SOURCE_DIR to override the default sibling checkout path." >&2 + echo "Set CHASTE_SOURCE_DIR to override the default checkout path." >&2 exit 1 fi @@ -27,5 +27,5 @@ mkdir -p "${output_dir}" "${build_dir}" export CHASTE_TEST_OUTPUT="${output_dir}" cd "${build_dir}" -nice -n 9 cmake "${chaste_dir}" +nice -n 9 cmake -DChaste_UPDATE_PROVENANCE=OFF "${chaste_dir}" diff --git a/chaste_sbml/SbmlRefModels/src/SbmlMath.hpp b/chaste_sbml/SbmlRefModels/src/SbmlMath.hpp index 35001701..51dc4494 100644 --- a/chaste_sbml/SbmlRefModels/src/SbmlMath.hpp +++ b/chaste_sbml/SbmlRefModels/src/SbmlMath.hpp @@ -95,6 +95,9 @@ constexpr bool or_(Args... args); // not_ inline bool not_(bool x) { return !x; } +// implies (a implies b == (not a) or b) +inline bool implies(double antecedent, double consequent) { return !antecedent || consequent; } + // xor_ template constexpr bool xor_(Args... args); @@ -149,17 +152,23 @@ inline double asech(double x) { return std::acosh(1.0 / x); } // factorial inline double factorial(double x) { return std::tgamma(x + 1.0); } -// max +// max (the max of a single value is that value) +constexpr double max(double only) { return only; } template constexpr double max(double first, double second, Args... rest); -// min +// min (the min of a single value is that value) +constexpr double min(double only) { return only; } template constexpr double min(double first, double second, Args... rest); // piecewise constexpr double piecewise(double otherwise); +// A piecewise with no otherwise clause: undefined (NaN) when the condition does not hold. Also the +// recursion base case for an even number of arguments (pieces without a trailing otherwise). +constexpr double piecewise(double value, bool condition); + constexpr double piecewise(double value, bool condition, double otherwise); template @@ -206,7 +215,9 @@ constexpr bool sbmlmath::or_(Args... args) template constexpr bool sbmlmath::xor_(Args... args) { - return (false ^ ... ^ args); + // Cast each operand to bool first: the arguments are numeric booleans (0/non-0), and ^ is + // bitwise (undefined for double), so xor them as bools -- true when an odd number are true. + return (false ^ ... ^ static_cast(args)); } // Relational (variadic templates) ============== @@ -296,6 +307,11 @@ constexpr double sbmlmath::piecewise(double otherwise) return otherwise; } +constexpr double sbmlmath::piecewise(double value, bool condition) +{ + return condition ? value : std::numeric_limits::quiet_NaN(); +} + constexpr double sbmlmath::piecewise(double value, bool condition, double otherwise) { return condition ? value : otherwise; diff --git a/chaste_sbml/SbmlRefModels/src/cases/README.md b/chaste_sbml/SbmlRefModels/src/cases/README.md index f23f57cb..951602bf 100644 --- a/chaste_sbml/SbmlRefModels/src/cases/README.md +++ b/chaste_sbml/SbmlRefModels/src/cases/README.md @@ -1,7 +1,10 @@ # SBML Test Suite Cases ## Semantic - -## Stochastic +Case support tracked in test/data/sbml_test_suite_status.csv ## Syntactic +Covered by libSBML + +## Stochastic +Not supported diff --git a/chaste_sbml/SbmlRefModels/src/reference/Chen2000/Chen2000SbmlOdeSystem.cpp b/chaste_sbml/SbmlRefModels/src/reference/Chen2000/Chen2000SbmlOdeSystem.cpp index ee1100fc..9681a281 100644 --- a/chaste_sbml/SbmlRefModels/src/reference/Chen2000/Chen2000SbmlOdeSystem.cpp +++ b/chaste_sbml/SbmlRefModels/src/reference/Chen2000/Chen2000SbmlOdeSystem.cpp @@ -64,119 +64,119 @@ void Chen2000SbmlOdeSystem::EvaluateYDerivatives(double time, const std::vector< void Chen2000SbmlOdeSystem::Initialise(double time) { - COMpartment = 1.0; // - Cln2 = 0.0078; // - ks_n2 = 0.0; // - ks_n2_ = 0.05; // - kd_n2 = 0.1; // - Clb2_T = 0.2342; // - Hct1_T = 1.0; // - ks_b2 = 0.002; // - ks_b2_ = 0.05; // - kd_b2 = 0.01; // - kd_b2_ = 2.0; // - kd_b2__ = 0.05; // - Clb5_T = 0.0614; // - ks_b5 = 0.006; // - ks_b5_ = 0.02; // - kd_b5 = 0.1; // - kd_b5_ = 0.25; // - Bck2_0 = 0.0027; // - Jn3 = 6.0; // - Dn3 = 1.0; // - Cln3_max = 0.02; // - Sic1_T = 0.1231; // - ks_c1 = 0.02; // - ks_c1_ = 0.1; // - Clb2_Sic1 = 0.079; // - kas_b2 = 50.0; // - kdi_b2 = 0.05; // - Clb5_Sic1 = 0.0207; // - kas_b5 = 50.0; // - kdi_b5 = 0.05; // - kd2_c1 = 0.3; // - epsilonc1_n3 = 20.0; // - epsilonc1_k2 = 2.0; // - epsilonc1_b5 = 1.0; // - epsilonc1_b2 = 0.067; // - Cdc20_T = 0.8332; // - ks_20 = 0.005; // - ks_20_ = 0.06; // - Cdc20 = 0.6848; // - ka_20 = 1.0; // - ki_20 = 0.1; // - ki_20_ = 10.0; // - Hct1 = 0.9946; // - ka_t1 = 0.04; // - ka_t1_ = 2.0; // - ki_t1 = 0.0; // - ki_t1_ = 0.64; // - Ji_t1 = 0.05; // - Ja_t1 = 0.05; // - epsiloni_t1_n2 = 1.0; // - epsiloni_t1_b5 = 0.5; // - epsiloni_t1_b2 = 1.0; // - mass = 0.6608; // - mu = 0.005776; // - ORI = 0.0; // - ks_ori = 2.0; // - kd_ori = 0.06; // - epsilonori_b2 = 0.4; // - BUD = 0.0; // - ks_bud = 0.3; // - kd_bud = 0.06; // - epsilonbud_b5 = 1.0; // - SPN = 0.0; // - ks_spn = 0.08; // - kd_spn = 0.06; // - J_spn = 0.2; // - ka_sbf = 1.0; // - ki_sbf = 0.5; // - ki_sbf_ = 6.0; // - Ji_sbf = 0.01; // - Ja_sbf = 0.01; // - epsilonsbf_n3 = 75.0; // - epsilonsbf_b5 = 0.5; // - ka_mcm = 1.0; // - ki_mcm = 0.15; // - Ji_mcm = 1.0; // - Ja_mcm = 1.0; // - ka_swi = 1.0; // - ki_swi = 0.3; // - ki_swi_ = 0.2; // - Ji_swi = 0.1; // - Ja_swi = 0.1; // - kd1_c1 = 0.01; // - kd_20 = 0.08; // - Jd2_c1 = 0.05; // - Vd_b2 = kd_b2 * (Hct1_T - Hct1) + kd_b2_ * Hct1 + kd_b2__ * Cdc20; // - Clb2 = Clb2_T - Clb2_Sic1; // - Clb5 = Clb5_T - Clb5_Sic1; // - Sic1 = Sic1_T - (Clb2_Sic1 + Clb5_Sic1); // - Vd_b5 = kd_b5 + kd_b5_ * Cdc20; // - Bck2 = Bck2_0 * mass; // - Cln3 = Cln3_max * Dn3 * mass / (Jn3 + Dn3 * mass); // - Vd2_c1 = kd2_c1 * (epsilonc1_n3 * Cln3 + epsilonc1_k2 * Bck2 + Cln2 + epsilonc1_b5 * Clb5 + epsilonc1_b2 * Clb2); // - Vi_20 = sm::piecewise(ki_20_, ORI >= 1.0, ki_20, SPN >= 1.0, 0.1); // - Vi_t1 = ki_t1 + ki_t1_ * (Cln3 + epsiloni_t1_n2 * Cln2 + epsiloni_t1_b5 * Clb5 + epsiloni_t1_b2 * Clb2); // - Va_sbf = ka_sbf * (Cln2 + epsilonsbf_n3 * (Cln3 + Bck2) + epsilonsbf_b5 * Clb5); // - SBF = 2.0 * Va_sbf * Ji_sbf / ((ki_sbf + ki_sbf_ * Clb2 + Va_sbf * Ji_sbf + (ki_sbf + ki_sbf_ * Clb2) * Ja_sbf - Va_sbf) + std::sqrt(std::pow((ki_sbf + ki_sbf_ * Clb2 + Va_sbf * Ji_sbf + (ki_sbf + ki_sbf_ * Clb2) * Ja_sbf - Va_sbf), 2.0) - 4.0 * Va_sbf * Ji_sbf * (ki_sbf + ki_sbf_ * Clb2 - Va_sbf))); // - MBF = SBF; // - Mcm1 = 2.0 * ka_mcm * Clb2 * Ji_mcm / ((ki_mcm + ka_mcm * Clb2 * Ji_mcm + ki_mcm * Ja_mcm - ka_mcm * Clb2) + std::sqrt(std::pow((ki_mcm + ka_mcm * Clb2 * Ji_mcm + ki_mcm * Ja_mcm - ka_mcm * Clb2), 2.0) - 4.0 * (ki_mcm - ka_mcm * Clb2) * ka_mcm * Clb2 * Ji_mcm)); // - Swi5 = 2.0 * ka_swi * Cdc20 * Ji_swi / ((ki_swi + ki_swi_ * Clb2 + ka_swi * Cdc20 * Ji_swi + (ki_swi + ki_swi_ * Clb2) * Ja_swi - ka_swi * Cdc20) + std::sqrt(std::pow((ki_swi + ki_swi_ * Clb2 + ka_swi * Cdc20 * Ji_swi + (ki_swi + ki_swi_ * Clb2) * Ja_swi - ka_swi * Cdc20), 2.0) - 4.0 * (ki_swi + ki_swi_ * Clb2 - ka_swi * Cdc20) * ka_swi * Cdc20 * Ji_swi)); // - d_Cln2_dt = mass * (ks_n2 + ks_n2_ * SBF) - kd_n2 * Cln2; // - d_Clb2_T_dt = mass * (ks_b2 + ks_b2_ * Mcm1) - Vd_b2 * Clb2_T; // - d_Clb5_T_dt = mass * (ks_b5 + ks_b5_ * MBF) - Vd_b5 * Clb5_T; // - d_Sic1_T_dt = ks_c1 + ks_c1_ * Swi5 - Sic1_T * (kd1_c1 + Vd2_c1 / (Jd2_c1 + Sic1_T)); // - d_Clb2_Sic1_dt = kas_b2 * Clb2 * Sic1 - Clb2_Sic1 * (kdi_b2 + Vd_b2 + kd1_c1 + Vd2_c1 / (Jd2_c1 + Sic1_T)); // - d_Clb5_Sic1_dt = kas_b5 * Clb5 * Sic1 - Clb5_Sic1 * (kdi_b5 + Vd_b5 + kd1_c1 + Vd2_c1 / (Jd2_c1 + Sic1_T)); // - d_Cdc20_T_dt = ks_20 + ks_20_ * Clb2 - kd_20 * Cdc20_T; // - d_Cdc20_dt = ka_20 * (Cdc20_T - Cdc20) - Cdc20 * (Vi_20 + kd_20); // - d_Hct1_dt = (ka_t1 + ka_t1_ * Cdc20) * (Hct1_T - Hct1) / (Ja_t1 + Hct1_T - Hct1) - Vi_t1 * Hct1 / (Ji_t1 + Hct1); // - d_mass_dt = mu * mass; // - d_ORI_dt = ks_ori * (Clb5 + epsilonori_b2 * Clb2) - kd_ori * ORI; // - d_BUD_dt = ks_bud * (Cln2 + Cln3 + epsilonbud_b5 * Clb5) - kd_bud * BUD; // - d_SPN_dt = ks_spn * Clb2 / (J_spn + Clb2) - kd_spn * SPN; // + COMpartment = 1.0; // + Cln2 = 0.0078; // + ks_n2 = 0.0; // + ks_n2_ = 0.05; // + kd_n2 = 0.1; // + Clb2_T = 0.2342; // + Hct1_T = 1.0; // + ks_b2 = 0.002; // + ks_b2_ = 0.05; // + kd_b2 = 0.01; // + kd_b2_ = 2.0; // + kd_b2__ = 0.05; // + Clb5_T = 0.0614; // + ks_b5 = 0.006; // + ks_b5_ = 0.02; // + kd_b5 = 0.1; // + kd_b5_ = 0.25; // + Bck2_0 = 0.0027; // + Jn3 = 6.0; // + Dn3 = 1.0; // + Cln3_max = 0.02; // + Sic1_T = 0.1231; // + ks_c1 = 0.02; // + ks_c1_ = 0.1; // + Clb2_Sic1 = 0.079; // + kas_b2 = 50.0; // + kdi_b2 = 0.05; // + Clb5_Sic1 = 0.0207; // + kas_b5 = 50.0; // + kdi_b5 = 0.05; // + kd2_c1 = 0.3; // + epsilonc1_n3 = 20.0; // + epsilonc1_k2 = 2.0; // + epsilonc1_b5 = 1.0; // + epsilonc1_b2 = 0.067; // + Cdc20_T = 0.8332; // + ks_20 = 0.005; // + ks_20_ = 0.06; // + Cdc20 = 0.6848; // + ka_20 = 1.0; // + ki_20 = 0.1; // + ki_20_ = 10.0; // + Hct1 = 0.9946; // + ka_t1 = 0.04; // + ka_t1_ = 2.0; // + ki_t1 = 0.0; // + ki_t1_ = 0.64; // + Ji_t1 = 0.05; // + Ja_t1 = 0.05; // + epsiloni_t1_n2 = 1.0; // + epsiloni_t1_b5 = 0.5; // + epsiloni_t1_b2 = 1.0; // + mass = 0.6608; // + mu = 0.005776; // + ORI = 0.0; // + ks_ori = 2.0; // + kd_ori = 0.06; // + epsilonori_b2 = 0.4; // + BUD = 0.0; // + ks_bud = 0.3; // + kd_bud = 0.06; // + epsilonbud_b5 = 1.0; // + SPN = 0.0; // + ks_spn = 0.08; // + kd_spn = 0.06; // + J_spn = 0.2; // + ka_sbf = 1.0; // + ki_sbf = 0.5; // + ki_sbf_ = 6.0; // + Ji_sbf = 0.01; // + Ja_sbf = 0.01; // + epsilonsbf_n3 = 75.0; // + epsilonsbf_b5 = 0.5; // + ka_mcm = 1.0; // + ki_mcm = 0.15; // + Ji_mcm = 1.0; // + Ja_mcm = 1.0; // + ka_swi = 1.0; // + ki_swi = 0.3; // + ki_swi_ = 0.2; // + Ji_swi = 0.1; // + Ja_swi = 0.1; // + kd1_c1 = 0.01; // + kd_20 = 0.08; // + Jd2_c1 = 0.05; // + Vd_b2 = kd_b2 * (Hct1_T - Hct1) + kd_b2_ * Hct1 + kd_b2__ * Cdc20; // + Clb2 = Clb2_T - Clb2_Sic1; // + Clb5 = Clb5_T - Clb5_Sic1; // + Sic1 = Sic1_T - (Clb2_Sic1 + Clb5_Sic1); // + Vd_b5 = kd_b5 + kd_b5_ * Cdc20; // + Bck2 = Bck2_0 * mass; // + Cln3 = Cln3_max * Dn3 * mass / (Jn3 + Dn3 * mass); // + Vd2_c1 = kd2_c1 * (epsilonc1_n3 * Cln3 + epsilonc1_k2 * Bck2 + Cln2 + epsilonc1_b5 * Clb5 + epsilonc1_b2 * Clb2); // + Vi_20 = sm::piecewise(ki_20_, ORI >= 1.0, ki_20, SPN >= 1.0, 0.1); // + Vi_t1 = ki_t1 + ki_t1_ * (Cln3 + epsiloni_t1_n2 * Cln2 + epsiloni_t1_b5 * Clb5 + epsiloni_t1_b2 * Clb2); // + Va_sbf = ka_sbf * (Cln2 + epsilonsbf_n3 * (Cln3 + Bck2) + epsilonsbf_b5 * Clb5); // + SBF = 2.0 * Va_sbf * Ji_sbf / ((ki_sbf + ki_sbf_ * Clb2 + Va_sbf * Ji_sbf + (ki_sbf + ki_sbf_ * Clb2) * Ja_sbf - Va_sbf) + std::sqrt(std::pow(ki_sbf + ki_sbf_ * Clb2 + Va_sbf * Ji_sbf + (ki_sbf + ki_sbf_ * Clb2) * Ja_sbf - Va_sbf, 2.0) - 4.0 * Va_sbf * Ji_sbf * (ki_sbf + ki_sbf_ * Clb2 - Va_sbf))); // + MBF = SBF; // + Mcm1 = 2.0 * ka_mcm * Clb2 * Ji_mcm / ((ki_mcm + ka_mcm * Clb2 * Ji_mcm + ki_mcm * Ja_mcm - ka_mcm * Clb2) + std::sqrt(std::pow(ki_mcm + ka_mcm * Clb2 * Ji_mcm + ki_mcm * Ja_mcm - ka_mcm * Clb2, 2.0) - 4.0 * (ki_mcm - ka_mcm * Clb2) * ka_mcm * Clb2 * Ji_mcm)); // + Swi5 = 2.0 * ka_swi * Cdc20 * Ji_swi / ((ki_swi + ki_swi_ * Clb2 + ka_swi * Cdc20 * Ji_swi + (ki_swi + ki_swi_ * Clb2) * Ja_swi - ka_swi * Cdc20) + std::sqrt(std::pow(ki_swi + ki_swi_ * Clb2 + ka_swi * Cdc20 * Ji_swi + (ki_swi + ki_swi_ * Clb2) * Ja_swi - ka_swi * Cdc20, 2.0) - 4.0 * (ki_swi + ki_swi_ * Clb2 - ka_swi * Cdc20) * ka_swi * Cdc20 * Ji_swi)); // + d_Cln2_dt = mass * (ks_n2 + ks_n2_ * SBF) - kd_n2 * Cln2; // + d_Clb2_T_dt = mass * (ks_b2 + ks_b2_ * Mcm1) - Vd_b2 * Clb2_T; // + d_Clb5_T_dt = mass * (ks_b5 + ks_b5_ * MBF) - Vd_b5 * Clb5_T; // + d_Sic1_T_dt = ks_c1 + ks_c1_ * Swi5 - Sic1_T * (kd1_c1 + Vd2_c1 / (Jd2_c1 + Sic1_T)); // + d_Clb2_Sic1_dt = kas_b2 * Clb2 * Sic1 - Clb2_Sic1 * (kdi_b2 + Vd_b2 + kd1_c1 + Vd2_c1 / (Jd2_c1 + Sic1_T)); // + d_Clb5_Sic1_dt = kas_b5 * Clb5 * Sic1 - Clb5_Sic1 * (kdi_b5 + Vd_b5 + kd1_c1 + Vd2_c1 / (Jd2_c1 + Sic1_T)); // + d_Cdc20_T_dt = ks_20 + ks_20_ * Clb2 - kd_20 * Cdc20_T; // + d_Cdc20_dt = ka_20 * (Cdc20_T - Cdc20) - Cdc20 * (Vi_20 + kd_20); // + d_Hct1_dt = (ka_t1 + ka_t1_ * Cdc20) * (Hct1_T - Hct1) / (Ja_t1 + Hct1_T - Hct1) - Vi_t1 * Hct1 / (Ji_t1 + Hct1); // + d_mass_dt = mu * mass; // + d_ORI_dt = ks_ori * (Clb5 + epsilonori_b2 * Clb2) - kd_ori * ORI; // + d_BUD_dt = ks_bud * (Cln2 + Cln3 + epsilonbud_b5 * Clb5) - kd_bud * BUD; // + d_SPN_dt = ks_spn * Clb2 / (J_spn + Clb2) - kd_spn * SPN; // mStateVariables.push_back(Cln2); mStateVariables.push_back(Clb2_T); @@ -392,34 +392,34 @@ std::vector Chen2000SbmlOdeSystem::RunModelEquations(double time, const kd_20 = GetParameter(69); Jd2_c1 = GetParameter(70); - Vd_b2 = kd_b2 * (Hct1_T - Hct1) + kd_b2_ * Hct1 + kd_b2__ * Cdc20; // - Clb2 = Clb2_T - Clb2_Sic1; // - Clb5 = Clb5_T - Clb5_Sic1; // - Sic1 = Sic1_T - (Clb2_Sic1 + Clb5_Sic1); // - Vd_b5 = kd_b5 + kd_b5_ * Cdc20; // - Bck2 = Bck2_0 * mass; // - Cln3 = Cln3_max * Dn3 * mass / (Jn3 + Dn3 * mass); // - Vd2_c1 = kd2_c1 * (epsilonc1_n3 * Cln3 + epsilonc1_k2 * Bck2 + Cln2 + epsilonc1_b5 * Clb5 + epsilonc1_b2 * Clb2); // - Vi_20 = sm::piecewise(ki_20_, ORI >= 1.0, ki_20, SPN >= 1.0, 0.1); // - Vi_t1 = ki_t1 + ki_t1_ * (Cln3 + epsiloni_t1_n2 * Cln2 + epsiloni_t1_b5 * Clb5 + epsiloni_t1_b2 * Clb2); // - Va_sbf = ka_sbf * (Cln2 + epsilonsbf_n3 * (Cln3 + Bck2) + epsilonsbf_b5 * Clb5); // - SBF = 2.0 * Va_sbf * Ji_sbf / ((ki_sbf + ki_sbf_ * Clb2 + Va_sbf * Ji_sbf + (ki_sbf + ki_sbf_ * Clb2) * Ja_sbf - Va_sbf) + std::sqrt(std::pow((ki_sbf + ki_sbf_ * Clb2 + Va_sbf * Ji_sbf + (ki_sbf + ki_sbf_ * Clb2) * Ja_sbf - Va_sbf), 2.0) - 4.0 * Va_sbf * Ji_sbf * (ki_sbf + ki_sbf_ * Clb2 - Va_sbf))); // - MBF = SBF; // - Mcm1 = 2.0 * ka_mcm * Clb2 * Ji_mcm / ((ki_mcm + ka_mcm * Clb2 * Ji_mcm + ki_mcm * Ja_mcm - ka_mcm * Clb2) + std::sqrt(std::pow((ki_mcm + ka_mcm * Clb2 * Ji_mcm + ki_mcm * Ja_mcm - ka_mcm * Clb2), 2.0) - 4.0 * (ki_mcm - ka_mcm * Clb2) * ka_mcm * Clb2 * Ji_mcm)); // - Swi5 = 2.0 * ka_swi * Cdc20 * Ji_swi / ((ki_swi + ki_swi_ * Clb2 + ka_swi * Cdc20 * Ji_swi + (ki_swi + ki_swi_ * Clb2) * Ja_swi - ka_swi * Cdc20) + std::sqrt(std::pow((ki_swi + ki_swi_ * Clb2 + ka_swi * Cdc20 * Ji_swi + (ki_swi + ki_swi_ * Clb2) * Ja_swi - ka_swi * Cdc20), 2.0) - 4.0 * (ki_swi + ki_swi_ * Clb2 - ka_swi * Cdc20) * ka_swi * Cdc20 * Ji_swi)); // - d_Cln2_dt = mass * (ks_n2 + ks_n2_ * SBF) - kd_n2 * Cln2; // - d_Clb2_T_dt = mass * (ks_b2 + ks_b2_ * Mcm1) - Vd_b2 * Clb2_T; // - d_Clb5_T_dt = mass * (ks_b5 + ks_b5_ * MBF) - Vd_b5 * Clb5_T; // - d_Sic1_T_dt = ks_c1 + ks_c1_ * Swi5 - Sic1_T * (kd1_c1 + Vd2_c1 / (Jd2_c1 + Sic1_T)); // - d_Clb2_Sic1_dt = kas_b2 * Clb2 * Sic1 - Clb2_Sic1 * (kdi_b2 + Vd_b2 + kd1_c1 + Vd2_c1 / (Jd2_c1 + Sic1_T)); // - d_Clb5_Sic1_dt = kas_b5 * Clb5 * Sic1 - Clb5_Sic1 * (kdi_b5 + Vd_b5 + kd1_c1 + Vd2_c1 / (Jd2_c1 + Sic1_T)); // - d_Cdc20_T_dt = ks_20 + ks_20_ * Clb2 - kd_20 * Cdc20_T; // - d_Cdc20_dt = ka_20 * (Cdc20_T - Cdc20) - Cdc20 * (Vi_20 + kd_20); // - d_Hct1_dt = (ka_t1 + ka_t1_ * Cdc20) * (Hct1_T - Hct1) / (Ja_t1 + Hct1_T - Hct1) - Vi_t1 * Hct1 / (Ji_t1 + Hct1); // - d_mass_dt = mu * mass; // - d_ORI_dt = ks_ori * (Clb5 + epsilonori_b2 * Clb2) - kd_ori * ORI; // - d_BUD_dt = ks_bud * (Cln2 + Cln3 + epsilonbud_b5 * Clb5) - kd_bud * BUD; // - d_SPN_dt = ks_spn * Clb2 / (J_spn + Clb2) - kd_spn * SPN; // + Vd_b2 = kd_b2 * (Hct1_T - Hct1) + kd_b2_ * Hct1 + kd_b2__ * Cdc20; // + Clb2 = Clb2_T - Clb2_Sic1; // + Clb5 = Clb5_T - Clb5_Sic1; // + Sic1 = Sic1_T - (Clb2_Sic1 + Clb5_Sic1); // + Vd_b5 = kd_b5 + kd_b5_ * Cdc20; // + Bck2 = Bck2_0 * mass; // + Cln3 = Cln3_max * Dn3 * mass / (Jn3 + Dn3 * mass); // + Vd2_c1 = kd2_c1 * (epsilonc1_n3 * Cln3 + epsilonc1_k2 * Bck2 + Cln2 + epsilonc1_b5 * Clb5 + epsilonc1_b2 * Clb2); // + Vi_20 = sm::piecewise(ki_20_, ORI >= 1.0, ki_20, SPN >= 1.0, 0.1); // + Vi_t1 = ki_t1 + ki_t1_ * (Cln3 + epsiloni_t1_n2 * Cln2 + epsiloni_t1_b5 * Clb5 + epsiloni_t1_b2 * Clb2); // + Va_sbf = ka_sbf * (Cln2 + epsilonsbf_n3 * (Cln3 + Bck2) + epsilonsbf_b5 * Clb5); // + SBF = 2.0 * Va_sbf * Ji_sbf / ((ki_sbf + ki_sbf_ * Clb2 + Va_sbf * Ji_sbf + (ki_sbf + ki_sbf_ * Clb2) * Ja_sbf - Va_sbf) + std::sqrt(std::pow(ki_sbf + ki_sbf_ * Clb2 + Va_sbf * Ji_sbf + (ki_sbf + ki_sbf_ * Clb2) * Ja_sbf - Va_sbf, 2.0) - 4.0 * Va_sbf * Ji_sbf * (ki_sbf + ki_sbf_ * Clb2 - Va_sbf))); // + MBF = SBF; // + Mcm1 = 2.0 * ka_mcm * Clb2 * Ji_mcm / ((ki_mcm + ka_mcm * Clb2 * Ji_mcm + ki_mcm * Ja_mcm - ka_mcm * Clb2) + std::sqrt(std::pow(ki_mcm + ka_mcm * Clb2 * Ji_mcm + ki_mcm * Ja_mcm - ka_mcm * Clb2, 2.0) - 4.0 * (ki_mcm - ka_mcm * Clb2) * ka_mcm * Clb2 * Ji_mcm)); // + Swi5 = 2.0 * ka_swi * Cdc20 * Ji_swi / ((ki_swi + ki_swi_ * Clb2 + ka_swi * Cdc20 * Ji_swi + (ki_swi + ki_swi_ * Clb2) * Ja_swi - ka_swi * Cdc20) + std::sqrt(std::pow(ki_swi + ki_swi_ * Clb2 + ka_swi * Cdc20 * Ji_swi + (ki_swi + ki_swi_ * Clb2) * Ja_swi - ka_swi * Cdc20, 2.0) - 4.0 * (ki_swi + ki_swi_ * Clb2 - ka_swi * Cdc20) * ka_swi * Cdc20 * Ji_swi)); // + d_Cln2_dt = mass * (ks_n2 + ks_n2_ * SBF) - kd_n2 * Cln2; // + d_Clb2_T_dt = mass * (ks_b2 + ks_b2_ * Mcm1) - Vd_b2 * Clb2_T; // + d_Clb5_T_dt = mass * (ks_b5 + ks_b5_ * MBF) - Vd_b5 * Clb5_T; // + d_Sic1_T_dt = ks_c1 + ks_c1_ * Swi5 - Sic1_T * (kd1_c1 + Vd2_c1 / (Jd2_c1 + Sic1_T)); // + d_Clb2_Sic1_dt = kas_b2 * Clb2 * Sic1 - Clb2_Sic1 * (kdi_b2 + Vd_b2 + kd1_c1 + Vd2_c1 / (Jd2_c1 + Sic1_T)); // + d_Clb5_Sic1_dt = kas_b5 * Clb5 * Sic1 - Clb5_Sic1 * (kdi_b5 + Vd_b5 + kd1_c1 + Vd2_c1 / (Jd2_c1 + Sic1_T)); // + d_Cdc20_T_dt = ks_20 + ks_20_ * Clb2 - kd_20 * Cdc20_T; // + d_Cdc20_dt = ka_20 * (Cdc20_T - Cdc20) - Cdc20 * (Vi_20 + kd_20); // + d_Hct1_dt = (ka_t1 + ka_t1_ * Cdc20) * (Hct1_T - Hct1) / (Ja_t1 + Hct1_T - Hct1) - Vi_t1 * Hct1 / (Ji_t1 + Hct1); // + d_mass_dt = mu * mass; // + d_ORI_dt = ks_ori * (Clb5 + epsilonori_b2 * Clb2) - kd_ori * ORI; // + d_BUD_dt = ks_bud * (Cln2 + Cln3 + epsilonbud_b5 * Clb5) - kd_bud * BUD; // + d_SPN_dt = ks_spn * Clb2 / (J_spn + Clb2) - kd_spn * SPN; // std::vector derivatives(13); derivatives[0] = d_Cln2_dt; diff --git a/chaste_sbml/SbmlRefModels/src/reference/Chen2004/Chen2004SbmlOdeSystem.cpp b/chaste_sbml/SbmlRefModels/src/reference/Chen2004/Chen2004SbmlOdeSystem.cpp index 9c1b0035..6355da1a 100644 --- a/chaste_sbml/SbmlRefModels/src/reference/Chen2004/Chen2004SbmlOdeSystem.cpp +++ b/chaste_sbml/SbmlRefModels/src/reference/Chen2004/Chen2004SbmlOdeSystem.cpp @@ -11,7 +11,7 @@ namespace sm = sbmlmath; Chen2004SbmlOdeSystem::Chen2004SbmlOdeSystem() - : AbstractSbmlOdeSystem(36, 143, 4) + : AbstractSbmlOdeSystem(36, 146, 4) { mpSystemInfo.reset(new CellwiseOdeSystemInformation); @@ -30,8 +30,8 @@ Chen2004SbmlOdeSystem::Chen2004SbmlOdeSystem() mEventSatisfied = { true, true, true, true }; // From SBML trigger initialValue mEventTriggered.resize(4, false); - mEventAdjustedParameters.resize(143, false); - mEventAdjustedParameterValues.resize(143, 0.0); + mEventAdjustedParameters.resize(146, false); + mEventAdjustedParameterValues.resize(146, 0.0); mEventAdjustedStateVars.resize(36, false); mEventAdjustedStateValues.resize(36, 0.0); @@ -104,7 +104,6 @@ std::vector Chen2004SbmlOdeSystem::ComputeDerivedQuantities(double time, dqs.push_back(cell); dqs.push_back(BCK2); - dqs.push_back(BUB2); dqs.push_back(CDC14T); dqs.push_back(CDC15i); dqs.push_back(CDC6T); @@ -113,8 +112,6 @@ std::vector Chen2004SbmlOdeSystem::ComputeDerivedQuantities(double time, dqs.push_back(CLB5T); dqs.push_back(CLN3); dqs.push_back(IE); - dqs.push_back(LTE1); - dqs.push_back(MAD2); dqs.push_back(MCM1); dqs.push_back(NET1T); dqs.push_back(PE); @@ -740,6 +737,9 @@ void Chen2004SbmlOdeSystem::Initialise(double time) SetDefaultInitialCondition(34, SWI5P); SetDefaultInitialCondition(35, TEM1GTP); + mParameters.push_back(BUB2); + mParameters.push_back(LTE1); + mParameters.push_back(MAD2); mParameters.push_back(b0); mParameters.push_back(bub2h); mParameters.push_back(bub2l); @@ -1076,8 +1076,23 @@ double Chen2004SbmlOdeSystem::ProcessModelEvents(double time, const std::vector< if (mEventTriggered[1]) { [[maybe_unused]] double event_priority = std::numeric_limits::max(); - MAD2 = mad2h; - BUB2 = bub2h; + // MAD2 = mad2h + if (!mEventAdjustedParameters[2] + || event_priority <= mEventAdjustedParameterPriority[2]) + { + mEventAdjustedParameters[2] = true; + mEventAdjustedParameterValues[2] = mad2h; + mEventAdjustedParameterPriority[2] = event_priority; + } + + // BUB2 = bub2h + if (!mEventAdjustedParameters[0] + || event_priority <= mEventAdjustedParameterPriority[0]) + { + mEventAdjustedParameters[0] = true; + mEventAdjustedParameterValues[0] = bub2h; + mEventAdjustedParameterPriority[0] = event_priority; + } } } @@ -1160,9 +1175,32 @@ double Chen2004SbmlOdeSystem::ProcessModelEvents(double time, const std::vector< if (mEventTriggered[2]) { [[maybe_unused]] double event_priority = std::numeric_limits::max(); - MAD2 = mad2l; - LTE1 = lte1h; - BUB2 = bub2l; + // MAD2 = mad2l + if (!mEventAdjustedParameters[2] + || event_priority <= mEventAdjustedParameterPriority[2]) + { + mEventAdjustedParameters[2] = true; + mEventAdjustedParameterValues[2] = mad2l; + mEventAdjustedParameterPriority[2] = event_priority; + } + + // LTE1 = lte1h + if (!mEventAdjustedParameters[1] + || event_priority <= mEventAdjustedParameterPriority[1]) + { + mEventAdjustedParameters[1] = true; + mEventAdjustedParameterValues[1] = lte1h; + mEventAdjustedParameterPriority[1] = event_priority; + } + + // BUB2 = bub2l + if (!mEventAdjustedParameters[0] + || event_priority <= mEventAdjustedParameterPriority[0]) + { + mEventAdjustedParameters[0] = true; + mEventAdjustedParameterValues[0] = bub2l; + mEventAdjustedParameterPriority[0] = event_priority; + } } } @@ -1254,7 +1292,15 @@ double Chen2004SbmlOdeSystem::ProcessModelEvents(double time, const std::vector< mEventAdjustedStatePriority[22] = event_priority; } - LTE1 = lte1l; + // LTE1 = lte1l + if (!mEventAdjustedParameters[1] + || event_priority <= mEventAdjustedParameterPriority[1]) + { + mEventAdjustedParameters[1] = true; + mEventAdjustedParameterValues[1] = lte1l; + mEventAdjustedParameterPriority[1] = event_priority; + } + // BUD = 0.0 if (!mEventAdjustedStateVars[0] || event_priority <= mEventAdjustedStatePriority[0]) @@ -1317,149 +1363,152 @@ std::vector Chen2004SbmlOdeSystem::RunModelEquations(double time, const SWI5P = rStateVariables[34]; TEM1GTP = rStateVariables[35]; - b0 = GetParameter(0); - bub2h = GetParameter(1); - bub2l = GetParameter(2); - C0 = GetParameter(3); - CDC15T = GetParameter(4); - Dn3 = GetParameter(5); - ebudb5 = GetParameter(6); - ebudn2 = GetParameter(7); - ebudn3 = GetParameter(8); - ec1b2 = GetParameter(9); - ec1b5 = GetParameter(10); - ec1k2 = GetParameter(11); - ec1n2 = GetParameter(12); - ec1n3 = GetParameter(13); - ef6b2 = GetParameter(14); - ef6b5 = GetParameter(15); - ef6k2 = GetParameter(16); - ef6n2 = GetParameter(17); - ef6n3 = GetParameter(18); - eicdhb2 = GetParameter(19); - eicdhb5 = GetParameter(20); - eicdhn2 = GetParameter(21); - eicdhn3 = GetParameter(22); - eorib2 = GetParameter(23); - eorib5 = GetParameter(24); - esbfb5 = GetParameter(25); - esbfn2 = GetParameter(26); - esbfn3 = GetParameter(27); - ESP1T = GetParameter(28); - IET = GetParameter(29); - J20ppx = GetParameter(30); - Jacdh = GetParameter(31); - Jaiep = GetParameter(32); - Jamcm = GetParameter(33); - Jasbf = GetParameter(34); - Jatem = GetParameter(35); - Jd2c1 = GetParameter(36); - Jd2f6 = GetParameter(37); - Jicdh = GetParameter(38); - Jiiep = GetParameter(39); - Jimcm = GetParameter(40); - Jisbf = GetParameter(41); - Jitem = GetParameter(42); - Jn3 = GetParameter(43); - Jpds = GetParameter(44); - Jspn = GetParameter(45); - ka15_p = GetParameter(46); - ka15_p_p = GetParameter(47); - ka15p = GetParameter(48); - ka20_p = GetParameter(49); - ka20_p_p = GetParameter(50); - kacdh_p = GetParameter(51); - kacdh_p_p = GetParameter(52); - kaiep = GetParameter(53); - kamcm = GetParameter(54); - kasb2 = GetParameter(55); - kasb5 = GetParameter(56); - kasbf = GetParameter(57); - kasesp = GetParameter(58); - kasf2 = GetParameter(59); - kasf5 = GetParameter(60); - kasrent = GetParameter(61); - kasrentp = GetParameter(62); - kaswi = GetParameter(63); - kd14 = GetParameter(64); - kd1c1 = GetParameter(65); - kd1f6 = GetParameter(66); - kd1pds_p = GetParameter(67); - kd20 = GetParameter(68); - kd2c1 = GetParameter(69); - kd2f6 = GetParameter(70); - kd2pds_p_p = GetParameter(71); - kd3c1 = GetParameter(72); - kd3f6 = GetParameter(73); - kd3pds_p_p = GetParameter(74); - kdb2_p = GetParameter(75); - kdb2_p_p = GetParameter(76); - kdb2p = GetParameter(77); - kdb5_p = GetParameter(78); - kdb5_p_p = GetParameter(79); - kdbud = GetParameter(80); - kdcdh = GetParameter(81); - kdib2 = GetParameter(82); - kdib5 = GetParameter(83); - kdiesp = GetParameter(84); - kdif2 = GetParameter(85); - kdif5 = GetParameter(86); - kdirent = GetParameter(87); - kdirentp = GetParameter(88); - kdn2 = GetParameter(89); - kdnet = GetParameter(90); - kdori = GetParameter(91); - kdppx_p = GetParameter(92); - kdppx_p_p = GetParameter(93); - kdspn = GetParameter(94); - kdswi = GetParameter(95); - KEZ = GetParameter(96); - KEZ2 = GetParameter(97); - ki15 = GetParameter(98); - kicdh_p = GetParameter(99); - kicdh_p_p = GetParameter(100); - kiiep = GetParameter(101); - kimcm = GetParameter(102); - kisbf_p = GetParameter(103); - kisbf_p_p = GetParameter(104); - kiswi = GetParameter(105); - kkpnet_p = GetParameter(106); - kkpnet_p_p = GetParameter(107); - kppc1 = GetParameter(108); - kppf6 = GetParameter(109); - kppnet_p = GetParameter(110); - kppnet_p_p = GetParameter(111); - ks14 = GetParameter(112); - ks1pds_p_p = GetParameter(113); - ks20_p = GetParameter(114); - ks20_p_p = GetParameter(115); - ks2pds_p_p = GetParameter(116); - ksb2_p = GetParameter(117); - ksb2_p_p = GetParameter(118); - ksb5_p = GetParameter(119); - ksb5_p_p = GetParameter(120); - ksbud = GetParameter(121); - ksc1_p = GetParameter(122); - ksc1_p_p = GetParameter(123); - kscdh = GetParameter(124); - ksf6_p = GetParameter(125); - ksf6_p_p = GetParameter(126); - ksf6_p_p_p = GetParameter(127); - ksn2_p = GetParameter(128); - ksn2_p_p = GetParameter(129); - ksnet = GetParameter(130); - ksori = GetParameter(131); - kspds_p = GetParameter(132); - ksppx = GetParameter(133); - ksspn = GetParameter(134); - ksswi_p = GetParameter(135); - ksswi_p_p = GetParameter(136); - lte1h = GetParameter(137); - lte1l = GetParameter(138); - mad2h = GetParameter(139); - mad2l = GetParameter(140); - mdt = GetParameter(141); - TEM1T = GetParameter(142); + BUB2 = GetParameter(0); + LTE1 = GetParameter(1); + MAD2 = GetParameter(2); + b0 = GetParameter(3); + bub2h = GetParameter(4); + bub2l = GetParameter(5); + C0 = GetParameter(6); + CDC15T = GetParameter(7); + Dn3 = GetParameter(8); + ebudb5 = GetParameter(9); + ebudn2 = GetParameter(10); + ebudn3 = GetParameter(11); + ec1b2 = GetParameter(12); + ec1b5 = GetParameter(13); + ec1k2 = GetParameter(14); + ec1n2 = GetParameter(15); + ec1n3 = GetParameter(16); + ef6b2 = GetParameter(17); + ef6b5 = GetParameter(18); + ef6k2 = GetParameter(19); + ef6n2 = GetParameter(20); + ef6n3 = GetParameter(21); + eicdhb2 = GetParameter(22); + eicdhb5 = GetParameter(23); + eicdhn2 = GetParameter(24); + eicdhn3 = GetParameter(25); + eorib2 = GetParameter(26); + eorib5 = GetParameter(27); + esbfb5 = GetParameter(28); + esbfn2 = GetParameter(29); + esbfn3 = GetParameter(30); + ESP1T = GetParameter(31); + IET = GetParameter(32); + J20ppx = GetParameter(33); + Jacdh = GetParameter(34); + Jaiep = GetParameter(35); + Jamcm = GetParameter(36); + Jasbf = GetParameter(37); + Jatem = GetParameter(38); + Jd2c1 = GetParameter(39); + Jd2f6 = GetParameter(40); + Jicdh = GetParameter(41); + Jiiep = GetParameter(42); + Jimcm = GetParameter(43); + Jisbf = GetParameter(44); + Jitem = GetParameter(45); + Jn3 = GetParameter(46); + Jpds = GetParameter(47); + Jspn = GetParameter(48); + ka15_p = GetParameter(49); + ka15_p_p = GetParameter(50); + ka15p = GetParameter(51); + ka20_p = GetParameter(52); + ka20_p_p = GetParameter(53); + kacdh_p = GetParameter(54); + kacdh_p_p = GetParameter(55); + kaiep = GetParameter(56); + kamcm = GetParameter(57); + kasb2 = GetParameter(58); + kasb5 = GetParameter(59); + kasbf = GetParameter(60); + kasesp = GetParameter(61); + kasf2 = GetParameter(62); + kasf5 = GetParameter(63); + kasrent = GetParameter(64); + kasrentp = GetParameter(65); + kaswi = GetParameter(66); + kd14 = GetParameter(67); + kd1c1 = GetParameter(68); + kd1f6 = GetParameter(69); + kd1pds_p = GetParameter(70); + kd20 = GetParameter(71); + kd2c1 = GetParameter(72); + kd2f6 = GetParameter(73); + kd2pds_p_p = GetParameter(74); + kd3c1 = GetParameter(75); + kd3f6 = GetParameter(76); + kd3pds_p_p = GetParameter(77); + kdb2_p = GetParameter(78); + kdb2_p_p = GetParameter(79); + kdb2p = GetParameter(80); + kdb5_p = GetParameter(81); + kdb5_p_p = GetParameter(82); + kdbud = GetParameter(83); + kdcdh = GetParameter(84); + kdib2 = GetParameter(85); + kdib5 = GetParameter(86); + kdiesp = GetParameter(87); + kdif2 = GetParameter(88); + kdif5 = GetParameter(89); + kdirent = GetParameter(90); + kdirentp = GetParameter(91); + kdn2 = GetParameter(92); + kdnet = GetParameter(93); + kdori = GetParameter(94); + kdppx_p = GetParameter(95); + kdppx_p_p = GetParameter(96); + kdspn = GetParameter(97); + kdswi = GetParameter(98); + KEZ = GetParameter(99); + KEZ2 = GetParameter(100); + ki15 = GetParameter(101); + kicdh_p = GetParameter(102); + kicdh_p_p = GetParameter(103); + kiiep = GetParameter(104); + kimcm = GetParameter(105); + kisbf_p = GetParameter(106); + kisbf_p_p = GetParameter(107); + kiswi = GetParameter(108); + kkpnet_p = GetParameter(109); + kkpnet_p_p = GetParameter(110); + kppc1 = GetParameter(111); + kppf6 = GetParameter(112); + kppnet_p = GetParameter(113); + kppnet_p_p = GetParameter(114); + ks14 = GetParameter(115); + ks1pds_p_p = GetParameter(116); + ks20_p = GetParameter(117); + ks20_p_p = GetParameter(118); + ks2pds_p_p = GetParameter(119); + ksb2_p = GetParameter(120); + ksb2_p_p = GetParameter(121); + ksb5_p = GetParameter(122); + ksb5_p_p = GetParameter(123); + ksbud = GetParameter(124); + ksc1_p = GetParameter(125); + ksc1_p_p = GetParameter(126); + kscdh = GetParameter(127); + ksf6_p = GetParameter(128); + ksf6_p_p = GetParameter(129); + ksf6_p_p_p = GetParameter(130); + ksn2_p = GetParameter(131); + ksn2_p_p = GetParameter(132); + ksnet = GetParameter(133); + ksori = GetParameter(134); + kspds_p = GetParameter(135); + ksppx = GetParameter(136); + ksspn = GetParameter(137); + ksswi_p = GetParameter(138); + ksswi_p_p = GetParameter(139); + lte1h = GetParameter(140); + lte1l = GetParameter(141); + mad2h = GetParameter(142); + mad2l = GetParameter(143); + mdt = GetParameter(144); + TEM1T = GetParameter(145); BCK2 = b0 * MASS; // CDC14T = CDC14 + RENT + RENTP; // @@ -1682,7 +1731,7 @@ inline double Chen2004SbmlOdeSystem::BB_218(double A1, double A2, double A3, dou inline double Chen2004SbmlOdeSystem::GK_219(double A1, double A2, double A3, double A4) { - return 2.0 * A4 * A1 / ((A2 - A1) + A3 * A2 + A4 * A1 + sm::root(2.0, std::pow(((A2 - A1) + A3 * A2 + A4 * A1), 2.0) - 4.0 * (A2 - A1) * A4 * A1)); + return 2.0 * A4 * A1 / ((A2 - A1) + A3 * A2 + A4 * A1 + sm::root(2.0, std::pow((A2 - A1) + A3 * A2 + A4 * A1, 2.0) - 4.0 * (A2 - A1) * A4 * A1)); } inline double Chen2004SbmlOdeSystem::MichaelisMenten_220(double M1, double J1, double k1, double S1) @@ -1855,9 +1904,6 @@ void CellwiseOdeSystemInformation::Initialise() this->mDerivedQuantityNames.push_back("BCK2"); this->mDerivedQuantityUnits.push_back("non-dim"); - this->mDerivedQuantityNames.push_back("BUB2"); - this->mDerivedQuantityUnits.push_back("non-dim"); - this->mDerivedQuantityNames.push_back("CDC14T"); this->mDerivedQuantityUnits.push_back("non-dim"); @@ -1882,12 +1928,6 @@ void CellwiseOdeSystemInformation::Initialise() this->mDerivedQuantityNames.push_back("IE"); this->mDerivedQuantityUnits.push_back("non-dim"); - this->mDerivedQuantityNames.push_back("LTE1"); - this->mDerivedQuantityUnits.push_back("non-dim"); - - this->mDerivedQuantityNames.push_back("MAD2"); - this->mDerivedQuantityUnits.push_back("non-dim"); - this->mDerivedQuantityNames.push_back("MCM1"); this->mDerivedQuantityUnits.push_back("non-dim"); @@ -2411,6 +2451,15 @@ void CellwiseOdeSystemInformation::Initialise() this->mDerivedQuantityUnits.push_back("non-dim"); // PARAMETERS + this->mParameterNames.push_back("BUB2"); + this->mParameterUnits.push_back("non-dim"); + + this->mParameterNames.push_back("LTE1"); + this->mParameterUnits.push_back("non-dim"); + + this->mParameterNames.push_back("MAD2"); + this->mParameterUnits.push_back("non-dim"); + this->mParameterNames.push_back("b0"); this->mParameterUnits.push_back("non-dim"); diff --git a/chaste_sbml/SbmlRefModels/src/reference/Chen2004/Chen2004SbmlOdeSystem.hpp b/chaste_sbml/SbmlRefModels/src/reference/Chen2004/Chen2004SbmlOdeSystem.hpp index 787accbe..47f083d9 100644 --- a/chaste_sbml/SbmlRefModels/src/reference/Chen2004/Chen2004SbmlOdeSystem.hpp +++ b/chaste_sbml/SbmlRefModels/src/reference/Chen2004/Chen2004SbmlOdeSystem.hpp @@ -27,6 +27,9 @@ class Chen2004SbmlOdeSystem : public AbstractSbmlOdeSystem } // PARAMETERS + double BUB2; // BUB2 + double LTE1; // LTE1 + double MAD2; // MAD2 double b0; // b0 double bub2h; // bub2h double bub2l; // bub2l @@ -249,7 +252,6 @@ class Chen2004SbmlOdeSystem : public AbstractSbmlOdeSystem // DERIVED QUANTITIES double cell; // cell double BCK2; // BCK2 - double BUB2; // BUB2 double CDC14T; // CDC14T double CDC15i; // CDC15i double CDC6T; // CDC6T @@ -258,8 +260,6 @@ class Chen2004SbmlOdeSystem : public AbstractSbmlOdeSystem double CLB5T; // CLB5T double CLN3; // CLN3 double IE; // IE - double LTE1; // LTE1 - double MAD2; // MAD2 double MCM1; // MCM1 double NET1T; // NET1T double PE; // PE diff --git a/chaste_sbml/SbmlRefModels/src/reference/Gardner1998/Gardner1998SbmlOdeSystem.cpp b/chaste_sbml/SbmlRefModels/src/reference/Gardner1998/Gardner1998SbmlOdeSystem.cpp index 2bee481e..78bf29b0 100644 --- a/chaste_sbml/SbmlRefModels/src/reference/Gardner1998/Gardner1998SbmlOdeSystem.cpp +++ b/chaste_sbml/SbmlRefModels/src/reference/Gardner1998/Gardner1998SbmlOdeSystem.cpp @@ -74,22 +74,22 @@ void Gardner1998SbmlOdeSystem::EvaluateYDerivatives(double time, const std::vect void Gardner1998SbmlOdeSystem::Initialise(double time) { - Cell = 1.0; // - C = 0.0; // - X = 0.0; // - M = 0.0; // - Y = 1.0; // - Z = 1.0; // - K6 = 0.3; // - V1p = 0.75; // - V3p = 0.3; // - C = C / Cell; // - X = X / Cell; // - M = M / Cell; // - Y = Y / Cell; // - Z = Z / Cell; // - V1 = C * V1p * std::pow((C + K6), -1.0); // - V3 = M * V3p; // + Cell = 1.0; // + C = 0.0; // + X = 0.0; // + M = 0.0; // + Y = 1.0; // + Z = 1.0; // + K6 = 0.3; // + V1p = 0.75; // + V3p = 0.3; // + C = C / Cell; // + X = X / Cell; // + M = M / Cell; // + Y = Y / Cell; // + Z = Z / Cell; // + V1 = C * V1p * std::pow(C + K6, -1.0); // + V3 = M * V3p; // // reaction1: { [[maybe_unused]] double vi = 0.1; @@ -105,7 +105,7 @@ void Gardner1998SbmlOdeSystem::Initialise(double time) // Qualify with this-> so the assignment targets the reaction member even when a local // parameter shadows its name (an SBML local parameter may shadow the reaction ID, e.g. a // reaction J1 with a local parameter also named J1). - this->reaction2 = C * k1 * X * std::pow((C + K5), -1.0); + this->reaction2 = C * k1 * X * std::pow(C + K5, -1.0); } // reaction3: { @@ -121,7 +121,7 @@ void Gardner1998SbmlOdeSystem::Initialise(double time) // Qualify with this-> so the assignment targets the reaction member even when a local // parameter shadows its name (an SBML local parameter may shadow the reaction ID, e.g. a // reaction J1 with a local parameter also named J1). - this->reaction4 = (1.0 + -1.0 * M) * V1 * std::pow((K1 + -1.0 * M + 1.0), -1.0); + this->reaction4 = (1.0 + -1.0 * M) * V1 * std::pow(K1 + -1.0 * M + 1.0, -1.0); } // reaction5: { @@ -130,7 +130,7 @@ void Gardner1998SbmlOdeSystem::Initialise(double time) // Qualify with this-> so the assignment targets the reaction member even when a local // parameter shadows its name (an SBML local parameter may shadow the reaction ID, e.g. a // reaction J1 with a local parameter also named J1). - this->reaction5 = M * V2 * std::pow((K2 + M), -1.0); + this->reaction5 = M * V2 * std::pow(K2 + M, -1.0); } // reaction6: { @@ -138,7 +138,7 @@ void Gardner1998SbmlOdeSystem::Initialise(double time) // Qualify with this-> so the assignment targets the reaction member even when a local // parameter shadows its name (an SBML local parameter may shadow the reaction ID, e.g. a // reaction J1 with a local parameter also named J1). - this->reaction6 = V3 * (1.0 + -1.0 * X) * std::pow((K3 + -1.0 * X + 1.0), -1.0); + this->reaction6 = V3 * (1.0 + -1.0 * X) * std::pow(K3 + -1.0 * X + 1.0, -1.0); } // reaction7: { @@ -147,7 +147,7 @@ void Gardner1998SbmlOdeSystem::Initialise(double time) // Qualify with this-> so the assignment targets the reaction member even when a local // parameter shadows its name (an SBML local parameter may shadow the reaction ID, e.g. a // reaction J1 with a local parameter also named J1). - this->reaction7 = V4 * X * std::pow((K4 + X), -1.0); + this->reaction7 = V4 * X * std::pow(K4 + X, -1.0); } // reaction8: { @@ -259,8 +259,8 @@ std::vector Gardner1998SbmlOdeSystem::RunModelEquations(double time, con V1p = GetParameter(1); V3p = GetParameter(2); - V1 = C * V1p * std::pow((C + K6), -1.0); // - V3 = M * V3p; // + V1 = C * V1p * std::pow(C + K6, -1.0); // + V3 = M * V3p; // // reaction1: { [[maybe_unused]] double vi = 0.1; @@ -276,7 +276,7 @@ std::vector Gardner1998SbmlOdeSystem::RunModelEquations(double time, con // Qualify with this-> so the assignment targets the reaction member even when a local // parameter shadows its name (an SBML local parameter may shadow the reaction ID, e.g. a // reaction J1 with a local parameter also named J1). - this->reaction2 = C * k1 * X * std::pow((C + K5), -1.0); + this->reaction2 = C * k1 * X * std::pow(C + K5, -1.0); } // reaction3: { @@ -292,7 +292,7 @@ std::vector Gardner1998SbmlOdeSystem::RunModelEquations(double time, con // Qualify with this-> so the assignment targets the reaction member even when a local // parameter shadows its name (an SBML local parameter may shadow the reaction ID, e.g. a // reaction J1 with a local parameter also named J1). - this->reaction4 = (1.0 + -1.0 * M) * V1 * std::pow((K1 + -1.0 * M + 1.0), -1.0); + this->reaction4 = (1.0 + -1.0 * M) * V1 * std::pow(K1 + -1.0 * M + 1.0, -1.0); } // reaction5: { @@ -301,7 +301,7 @@ std::vector Gardner1998SbmlOdeSystem::RunModelEquations(double time, con // Qualify with this-> so the assignment targets the reaction member even when a local // parameter shadows its name (an SBML local parameter may shadow the reaction ID, e.g. a // reaction J1 with a local parameter also named J1). - this->reaction5 = M * V2 * std::pow((K2 + M), -1.0); + this->reaction5 = M * V2 * std::pow(K2 + M, -1.0); } // reaction6: { @@ -309,7 +309,7 @@ std::vector Gardner1998SbmlOdeSystem::RunModelEquations(double time, con // Qualify with this-> so the assignment targets the reaction member even when a local // parameter shadows its name (an SBML local parameter may shadow the reaction ID, e.g. a // reaction J1 with a local parameter also named J1). - this->reaction6 = V3 * (1.0 + -1.0 * X) * std::pow((K3 + -1.0 * X + 1.0), -1.0); + this->reaction6 = V3 * (1.0 + -1.0 * X) * std::pow(K3 + -1.0 * X + 1.0, -1.0); } // reaction7: { @@ -318,7 +318,7 @@ std::vector Gardner1998SbmlOdeSystem::RunModelEquations(double time, con // Qualify with this-> so the assignment targets the reaction member even when a local // parameter shadows its name (an SBML local parameter may shadow the reaction ID, e.g. a // reaction J1 with a local parameter also named J1). - this->reaction7 = V4 * X * std::pow((K4 + X), -1.0); + this->reaction7 = V4 * X * std::pow(K4 + X, -1.0); } // reaction8: { diff --git a/chaste_sbml/SbmlRefModels/src/reference/Goldbeter1991/Goldbeter1991SbmlOdeSystem.cpp b/chaste_sbml/SbmlRefModels/src/reference/Goldbeter1991/Goldbeter1991SbmlOdeSystem.cpp index 4cc43c96..65deeb9a 100644 --- a/chaste_sbml/SbmlRefModels/src/reference/Goldbeter1991/Goldbeter1991SbmlOdeSystem.cpp +++ b/chaste_sbml/SbmlRefModels/src/reference/Goldbeter1991/Goldbeter1991SbmlOdeSystem.cpp @@ -64,15 +64,15 @@ void Goldbeter1991SbmlOdeSystem::EvaluateYDerivatives(double time, const std::ve void Goldbeter1991SbmlOdeSystem::Initialise(double time) { - cell = 1.0; // - C = 0.01; // - M = 0.01; // - X = 0.01; // - VM1 = 3.0; // - VM3 = 1.0; // - Kc = 0.5; // - V1 = C * VM1 * std::pow((C + Kc), -1.0); // - V3 = M * VM3; // + cell = 1.0; // + C = 0.01; // + M = 0.01; // + X = 0.01; // + VM1 = 3.0; // + VM3 = 1.0; // + Kc = 0.5; // + V1 = C * VM1 * std::pow(C + Kc, -1.0); // + V3 = M * VM3; // // reaction1: { [[maybe_unused]] double vi = 0.025; @@ -96,7 +96,7 @@ void Goldbeter1991SbmlOdeSystem::Initialise(double time) // Qualify with this-> so the assignment targets the reaction member even when a local // parameter shadows its name (an SBML local parameter may shadow the reaction ID, e.g. a // reaction J1 with a local parameter also named J1). - this->reaction3 = C * cell * vd * X * std::pow((C + Kd), -1.0); + this->reaction3 = C * cell * vd * X * std::pow(C + Kd, -1.0); } // reaction4: { @@ -104,7 +104,7 @@ void Goldbeter1991SbmlOdeSystem::Initialise(double time) // Qualify with this-> so the assignment targets the reaction member even when a local // parameter shadows its name (an SBML local parameter may shadow the reaction ID, e.g. a // reaction J1 with a local parameter also named J1). - this->reaction4 = cell * (1.0 + -1.0 * M) * V1 * std::pow((K1 + -1.0 * M + 1.0), -1.0); + this->reaction4 = cell * (1.0 + -1.0 * M) * V1 * std::pow(K1 + -1.0 * M + 1.0, -1.0); } // reaction5: { @@ -113,7 +113,7 @@ void Goldbeter1991SbmlOdeSystem::Initialise(double time) // Qualify with this-> so the assignment targets the reaction member even when a local // parameter shadows its name (an SBML local parameter may shadow the reaction ID, e.g. a // reaction J1 with a local parameter also named J1). - this->reaction5 = cell * M * V2 * std::pow((K2 + M), -1.0); + this->reaction5 = cell * M * V2 * std::pow(K2 + M, -1.0); } // reaction6: { @@ -121,7 +121,7 @@ void Goldbeter1991SbmlOdeSystem::Initialise(double time) // Qualify with this-> so the assignment targets the reaction member even when a local // parameter shadows its name (an SBML local parameter may shadow the reaction ID, e.g. a // reaction J1 with a local parameter also named J1). - this->reaction6 = cell * V3 * (1.0 + -1.0 * X) * std::pow((K3 + -1.0 * X + 1.0), -1.0); + this->reaction6 = cell * V3 * (1.0 + -1.0 * X) * std::pow(K3 + -1.0 * X + 1.0, -1.0); } // reaction7: { @@ -130,7 +130,7 @@ void Goldbeter1991SbmlOdeSystem::Initialise(double time) // Qualify with this-> so the assignment targets the reaction member even when a local // parameter shadows its name (an SBML local parameter may shadow the reaction ID, e.g. a // reaction J1 with a local parameter also named J1). - this->reaction7 = cell * V4 * X * std::pow((K4 + X), -1.0); + this->reaction7 = cell * V4 * X * std::pow(K4 + X, -1.0); } d_C_dt = (reaction1 - reaction2 - reaction3) / cell; // d_M_dt = (reaction4 - reaction5) / cell; // @@ -184,8 +184,8 @@ std::vector Goldbeter1991SbmlOdeSystem::RunModelEquations(double time, c VM3 = GetParameter(1); Kc = GetParameter(2); - V1 = C * VM1 * std::pow((C + Kc), -1.0); // - V3 = M * VM3; // + V1 = C * VM1 * std::pow(C + Kc, -1.0); // + V3 = M * VM3; // // reaction1: { [[maybe_unused]] double vi = 0.025; @@ -209,7 +209,7 @@ std::vector Goldbeter1991SbmlOdeSystem::RunModelEquations(double time, c // Qualify with this-> so the assignment targets the reaction member even when a local // parameter shadows its name (an SBML local parameter may shadow the reaction ID, e.g. a // reaction J1 with a local parameter also named J1). - this->reaction3 = C * cell * vd * X * std::pow((C + Kd), -1.0); + this->reaction3 = C * cell * vd * X * std::pow(C + Kd, -1.0); } // reaction4: { @@ -217,7 +217,7 @@ std::vector Goldbeter1991SbmlOdeSystem::RunModelEquations(double time, c // Qualify with this-> so the assignment targets the reaction member even when a local // parameter shadows its name (an SBML local parameter may shadow the reaction ID, e.g. a // reaction J1 with a local parameter also named J1). - this->reaction4 = cell * (1.0 + -1.0 * M) * V1 * std::pow((K1 + -1.0 * M + 1.0), -1.0); + this->reaction4 = cell * (1.0 + -1.0 * M) * V1 * std::pow(K1 + -1.0 * M + 1.0, -1.0); } // reaction5: { @@ -226,7 +226,7 @@ std::vector Goldbeter1991SbmlOdeSystem::RunModelEquations(double time, c // Qualify with this-> so the assignment targets the reaction member even when a local // parameter shadows its name (an SBML local parameter may shadow the reaction ID, e.g. a // reaction J1 with a local parameter also named J1). - this->reaction5 = cell * M * V2 * std::pow((K2 + M), -1.0); + this->reaction5 = cell * M * V2 * std::pow(K2 + M, -1.0); } // reaction6: { @@ -234,7 +234,7 @@ std::vector Goldbeter1991SbmlOdeSystem::RunModelEquations(double time, c // Qualify with this-> so the assignment targets the reaction member even when a local // parameter shadows its name (an SBML local parameter may shadow the reaction ID, e.g. a // reaction J1 with a local parameter also named J1). - this->reaction6 = cell * V3 * (1.0 + -1.0 * X) * std::pow((K3 + -1.0 * X + 1.0), -1.0); + this->reaction6 = cell * V3 * (1.0 + -1.0 * X) * std::pow(K3 + -1.0 * X + 1.0, -1.0); } // reaction7: { @@ -243,7 +243,7 @@ std::vector Goldbeter1991SbmlOdeSystem::RunModelEquations(double time, c // Qualify with this-> so the assignment targets the reaction member even when a local // parameter shadows its name (an SBML local parameter may shadow the reaction ID, e.g. a // reaction J1 with a local parameter also named J1). - this->reaction7 = cell * V4 * X * std::pow((K4 + X), -1.0); + this->reaction7 = cell * V4 * X * std::pow(K4 + X, -1.0); } d_C_dt = (reaction1 - reaction2 - reaction3) / cell; // d_M_dt = (reaction4 - reaction5) / cell; // diff --git a/chaste_sbml/SbmlRefModels/src/reference/TysonNovak2001/TysonNovak2001SbmlOdeSystem.cpp b/chaste_sbml/SbmlRefModels/src/reference/TysonNovak2001/TysonNovak2001SbmlOdeSystem.cpp index ba1391b3..cb94a31f 100644 --- a/chaste_sbml/SbmlRefModels/src/reference/TysonNovak2001/TysonNovak2001SbmlOdeSystem.cpp +++ b/chaste_sbml/SbmlRefModels/src/reference/TysonNovak2001/TysonNovak2001SbmlOdeSystem.cpp @@ -109,83 +109,83 @@ void TysonNovak2001SbmlOdeSystem::EvaluateYDerivatives(double time, const std::v void TysonNovak2001SbmlOdeSystem::Initialise(double time) { - cell = 1.0; // - CycBt = 0.001; // - Cdc20a = 0.001; // - Cdh1 = 0.001; // - m = 0.5; // - Cdc20t = 0.001; // - IEP = 0.001; // - CKIt = 0.001; // - SK = 0.001; // - k1 = 0.04; // - k2p = 0.04; // - k2pp = 1.0; // - k2ppp = 1.0; // - k3p = 1.0; // - k3pp = 10.0; // - J3 = 0.04; // - k4 = 35.0; // - k5p = 0.005; // - k5pp = 0.2; // - J5 = 0.3; // - k6 = 0.1; // - n = 4.0; // - k7 = 1.0; // - J7 = 0.001; // - k8 = 0.5; // - J8 = 0.001; // - k9 = 0.1; // - k10 = 0.02; // - mu = 0.005; // - k11 = 1.0; // - k12p = 0.2; // - k12pp = 50.0; // - mmax = 10.0; // - k12ppp = 100.0; // - Keq = 1000.0; // - k13 = 1.0; // - k14 = 1.0; // - k15p = 1.5; // - k15pp = 0.05; // - k16p = 1.0; // - k16pp = 3.0; // - J15 = 0.01; // - J16 = 0.01; // - k4p = 2.0; // - J4 = 0.04; // - CycB = CycBt - 2.0 * CycBt * CKIt / (CycBt + CKIt + 1.0 / Keq + std::pow((std::pow((CycBt + CKIt + 1.0 / Keq), 2.0) - 4.0 * CycBt * CKIt), (1.0 / 2.0))); // - Trimer = 2.0 * CycBt * CKIt / (CycBt + CKIt + 1.0 / Keq + std::pow((std::pow((CycBt + CKIt + 1.0 / Keq), 2.0) - 4.0 * CycBt * CKIt), (1.0 / 2.0))); // - Mad = 1.0; // - TF = GK(k15p * m + k15pp * SK, k16p + k16pp * m * CycB, J15, J16); // - CycBt_synthesis = k1; // - CycBdegradation = k2p * CycBt; // - CycBdegradationviaCdh1 = k2pp * Cdh1 * CycBt; // - CycBtdegradationviaCdc20a = k2ppp * Cdc20a * CycBt; // - Cdh1synthesis = (k3p + k3pp * Cdc20a) * (1.0 - Cdh1) / (J3 + 1.0 - Cdh1); // - Cdh1degradation = (k4p * SK * Cdh1 + k4 * m * CycB * Cdh1) / (J4 + Cdh1); // - Cdc20tsynthesis = k5p + k5pp * std::pow((CycB * m / J5), n) / (1.0 + std::pow((CycB * m / J5), n)); // - Cdc20t_deg = k6 * Cdc20t; // - Cdc20activation = k7 * IEP * (Cdc20t - Cdc20a) / (J7 + Cdc20t - Cdc20a); // - Cdc20ainhibition = k8 * Mad * Cdc20a / (J8 + Cdc20a); // - Cdc20adegradation = k6 * Cdc20a; // - IEPsynthesis = k9 * m * CycB * (1.0 - IEP); // - IEPdegradation = k10 * IEP; // - growth = mu * m * (1.0 - m / mmax); // - CKItsynthesis = k11; // - CKIdegradation = k12p * CKIt; // - CKItphosphorilationviaSK = k12pp * SK * CKIt; // - eq_7 = k12ppp * m * CycB * CKIt; // - SKsynthesis = k13 * TF; // - SKdegradation = k14 * SK; // - d_CycBt_dt = CycBt_synthesis - CycBdegradation - CycBdegradationviaCdh1 - CycBtdegradationviaCdc20a; // - d_Cdc20a_dt = Cdc20activation - Cdc20ainhibition - Cdc20adegradation; // - d_Cdh1_dt = Cdh1synthesis - Cdh1degradation; // - d_m_dt = growth; // - d_Cdc20t_dt = Cdc20tsynthesis - Cdc20t_deg; // - d_IEP_dt = IEPsynthesis - IEPdegradation; // - d_CKIt_dt = CKItsynthesis - CKIdegradation - CKItphosphorilationviaSK - eq_7; // - d_SK_dt = SKsynthesis - SKdegradation; // + cell = 1.0; // + CycBt = 0.001; // + Cdc20a = 0.001; // + Cdh1 = 0.001; // + m = 0.5; // + Cdc20t = 0.001; // + IEP = 0.001; // + CKIt = 0.001; // + SK = 0.001; // + k1 = 0.04; // + k2p = 0.04; // + k2pp = 1.0; // + k2ppp = 1.0; // + k3p = 1.0; // + k3pp = 10.0; // + J3 = 0.04; // + k4 = 35.0; // + k5p = 0.005; // + k5pp = 0.2; // + J5 = 0.3; // + k6 = 0.1; // + n = 4.0; // + k7 = 1.0; // + J7 = 0.001; // + k8 = 0.5; // + J8 = 0.001; // + k9 = 0.1; // + k10 = 0.02; // + mu = 0.005; // + k11 = 1.0; // + k12p = 0.2; // + k12pp = 50.0; // + mmax = 10.0; // + k12ppp = 100.0; // + Keq = 1000.0; // + k13 = 1.0; // + k14 = 1.0; // + k15p = 1.5; // + k15pp = 0.05; // + k16p = 1.0; // + k16pp = 3.0; // + J15 = 0.01; // + J16 = 0.01; // + k4p = 2.0; // + J4 = 0.04; // + CycB = CycBt - 2.0 * CycBt * CKIt / (CycBt + CKIt + 1.0 / Keq + std::pow(std::pow(CycBt + CKIt + 1.0 / Keq, 2.0) - 4.0 * CycBt * CKIt, 1.0 / 2.0)); // + Trimer = 2.0 * CycBt * CKIt / (CycBt + CKIt + 1.0 / Keq + std::pow(std::pow(CycBt + CKIt + 1.0 / Keq, 2.0) - 4.0 * CycBt * CKIt, 1.0 / 2.0)); // + Mad = 1.0; // + TF = GK(k15p * m + k15pp * SK, k16p + k16pp * m * CycB, J15, J16); // + CycBt_synthesis = k1; // + CycBdegradation = k2p * CycBt; // + CycBdegradationviaCdh1 = k2pp * Cdh1 * CycBt; // + CycBtdegradationviaCdc20a = k2ppp * Cdc20a * CycBt; // + Cdh1synthesis = (k3p + k3pp * Cdc20a) * (1.0 - Cdh1) / (J3 + 1.0 - Cdh1); // + Cdh1degradation = (k4p * SK * Cdh1 + k4 * m * CycB * Cdh1) / (J4 + Cdh1); // + Cdc20tsynthesis = k5p + k5pp * std::pow(CycB * m / J5, n) / (1.0 + std::pow(CycB * m / J5, n)); // + Cdc20t_deg = k6 * Cdc20t; // + Cdc20activation = k7 * IEP * (Cdc20t - Cdc20a) / (J7 + Cdc20t - Cdc20a); // + Cdc20ainhibition = k8 * Mad * Cdc20a / (J8 + Cdc20a); // + Cdc20adegradation = k6 * Cdc20a; // + IEPsynthesis = k9 * m * CycB * (1.0 - IEP); // + IEPdegradation = k10 * IEP; // + growth = mu * m * (1.0 - m / mmax); // + CKItsynthesis = k11; // + CKIdegradation = k12p * CKIt; // + CKItphosphorilationviaSK = k12pp * SK * CKIt; // + eq_7 = k12ppp * m * CycB * CKIt; // + SKsynthesis = k13 * TF; // + SKdegradation = k14 * SK; // + d_CycBt_dt = CycBt_synthesis - CycBdegradation - CycBdegradationviaCdh1 - CycBtdegradationviaCdc20a; // + d_Cdc20a_dt = Cdc20activation - Cdc20ainhibition - Cdc20adegradation; // + d_Cdh1_dt = Cdh1synthesis - Cdh1degradation; // + d_m_dt = growth; // + d_Cdc20t_dt = Cdc20tsynthesis - Cdc20t_deg; // + d_IEP_dt = IEPsynthesis - IEPdegradation; // + d_CKIt_dt = CKItsynthesis - CKIdegradation - CKItphosphorilationviaSK - eq_7; // + d_SK_dt = SKsynthesis - SKdegradation; // mStateVariables.push_back(CycBt); mStateVariables.push_back(Cdc20a); @@ -406,38 +406,38 @@ std::vector TysonNovak2001SbmlOdeSystem::RunModelEquations(double time, k4p = GetParameter(34); J4 = GetParameter(35); - CycB = CycBt - 2.0 * CycBt * CKIt / (CycBt + CKIt + 1.0 / Keq + std::pow((std::pow((CycBt + CKIt + 1.0 / Keq), 2.0) - 4.0 * CycBt * CKIt), (1.0 / 2.0))); // - Trimer = 2.0 * CycBt * CKIt / (CycBt + CKIt + 1.0 / Keq + std::pow((std::pow((CycBt + CKIt + 1.0 / Keq), 2.0) - 4.0 * CycBt * CKIt), (1.0 / 2.0))); // - Mad = 1.0; // - TF = GK(k15p * m + k15pp * SK, k16p + k16pp * m * CycB, J15, J16); // - CycBt_synthesis = k1; // - CycBdegradation = k2p * CycBt; // - CycBdegradationviaCdh1 = k2pp * Cdh1 * CycBt; // - CycBtdegradationviaCdc20a = k2ppp * Cdc20a * CycBt; // - Cdh1synthesis = (k3p + k3pp * Cdc20a) * (1.0 - Cdh1) / (J3 + 1.0 - Cdh1); // - Cdh1degradation = (k4p * SK * Cdh1 + k4 * m * CycB * Cdh1) / (J4 + Cdh1); // - Cdc20tsynthesis = k5p + k5pp * std::pow((CycB * m / J5), n) / (1.0 + std::pow((CycB * m / J5), n)); // - Cdc20t_deg = k6 * Cdc20t; // - Cdc20activation = k7 * IEP * (Cdc20t - Cdc20a) / (J7 + Cdc20t - Cdc20a); // - Cdc20ainhibition = k8 * Mad * Cdc20a / (J8 + Cdc20a); // - Cdc20adegradation = k6 * Cdc20a; // - IEPsynthesis = k9 * m * CycB * (1.0 - IEP); // - IEPdegradation = k10 * IEP; // - growth = mu * m * (1.0 - m / mmax); // - CKItsynthesis = k11; // - CKIdegradation = k12p * CKIt; // - CKItphosphorilationviaSK = k12pp * SK * CKIt; // - eq_7 = k12ppp * m * CycB * CKIt; // - SKsynthesis = k13 * TF; // - SKdegradation = k14 * SK; // - d_CycBt_dt = CycBt_synthesis - CycBdegradation - CycBdegradationviaCdh1 - CycBtdegradationviaCdc20a; // - d_Cdc20a_dt = Cdc20activation - Cdc20ainhibition - Cdc20adegradation; // - d_Cdh1_dt = Cdh1synthesis - Cdh1degradation; // - d_m_dt = growth; // - d_Cdc20t_dt = Cdc20tsynthesis - Cdc20t_deg; // - d_IEP_dt = IEPsynthesis - IEPdegradation; // - d_CKIt_dt = CKItsynthesis - CKIdegradation - CKItphosphorilationviaSK - eq_7; // - d_SK_dt = SKsynthesis - SKdegradation; // + CycB = CycBt - 2.0 * CycBt * CKIt / (CycBt + CKIt + 1.0 / Keq + std::pow(std::pow(CycBt + CKIt + 1.0 / Keq, 2.0) - 4.0 * CycBt * CKIt, 1.0 / 2.0)); // + Trimer = 2.0 * CycBt * CKIt / (CycBt + CKIt + 1.0 / Keq + std::pow(std::pow(CycBt + CKIt + 1.0 / Keq, 2.0) - 4.0 * CycBt * CKIt, 1.0 / 2.0)); // + Mad = 1.0; // + TF = GK(k15p * m + k15pp * SK, k16p + k16pp * m * CycB, J15, J16); // + CycBt_synthesis = k1; // + CycBdegradation = k2p * CycBt; // + CycBdegradationviaCdh1 = k2pp * Cdh1 * CycBt; // + CycBtdegradationviaCdc20a = k2ppp * Cdc20a * CycBt; // + Cdh1synthesis = (k3p + k3pp * Cdc20a) * (1.0 - Cdh1) / (J3 + 1.0 - Cdh1); // + Cdh1degradation = (k4p * SK * Cdh1 + k4 * m * CycB * Cdh1) / (J4 + Cdh1); // + Cdc20tsynthesis = k5p + k5pp * std::pow(CycB * m / J5, n) / (1.0 + std::pow(CycB * m / J5, n)); // + Cdc20t_deg = k6 * Cdc20t; // + Cdc20activation = k7 * IEP * (Cdc20t - Cdc20a) / (J7 + Cdc20t - Cdc20a); // + Cdc20ainhibition = k8 * Mad * Cdc20a / (J8 + Cdc20a); // + Cdc20adegradation = k6 * Cdc20a; // + IEPsynthesis = k9 * m * CycB * (1.0 - IEP); // + IEPdegradation = k10 * IEP; // + growth = mu * m * (1.0 - m / mmax); // + CKItsynthesis = k11; // + CKIdegradation = k12p * CKIt; // + CKItphosphorilationviaSK = k12pp * SK * CKIt; // + eq_7 = k12ppp * m * CycB * CKIt; // + SKsynthesis = k13 * TF; // + SKdegradation = k14 * SK; // + d_CycBt_dt = CycBt_synthesis - CycBdegradation - CycBdegradationviaCdh1 - CycBtdegradationviaCdc20a; // + d_Cdc20a_dt = Cdc20activation - Cdc20ainhibition - Cdc20adegradation; // + d_Cdh1_dt = Cdh1synthesis - Cdh1degradation; // + d_m_dt = growth; // + d_Cdc20t_dt = Cdc20tsynthesis - Cdc20t_deg; // + d_IEP_dt = IEPsynthesis - IEPdegradation; // + d_CKIt_dt = CKItsynthesis - CKIdegradation - CKItphosphorilationviaSK - eq_7; // + d_SK_dt = SKsynthesis - SKdegradation; // std::vector derivatives(8); derivatives[0] = d_CycBt_dt; @@ -454,7 +454,7 @@ std::vector TysonNovak2001SbmlOdeSystem::RunModelEquations(double time, // MODEL FUNCTIONS inline double TysonNovak2001SbmlOdeSystem::GK(double A1, double A2, double A3, double A4) { - return 2.0 * A4 * A1 / ((A2 - A1) + A3 * A2 + A4 * A1 + sm::root(2.0, std::pow(((A2 - A1) + A3 * A2 + A4 * A1), 2.0) - 4.0 * (A2 - A1) * A4 * A1)); + return 2.0 * A4 * A1 / ((A2 - A1) + A3 * A2 + A4 * A1 + sm::root(2.0, std::pow((A2 - A1) + A3 * A2 + A4 * A1, 2.0) - 4.0 * (A2 - A1) * A4 * A1)); } template <> diff --git a/chaste_sbml/SbmlRefModels/test/CMakeLists.txt b/chaste_sbml/SbmlRefModels/test/CMakeLists.txt index cea25f65..abee9aea 100644 --- a/chaste_sbml/SbmlRefModels/test/CMakeLists.txt +++ b/chaste_sbml/SbmlRefModels/test/CMakeLists.txt @@ -1 +1,15 @@ chaste_do_test_project(SbmlRefModels) + +# Precompile the heavy Chaste headers shared by every generated case test (the CVODE solver and +# test helpers), to speed up rebuilds. FakePetscSetup and the per-test model header are left out. +get_property(_sbml_test_targets DIRECTORY PROPERTY BUILDSYSTEM_TARGETS) +foreach(_sbml_test_target ${_sbml_test_targets}) + get_target_property(_sbml_test_type ${_sbml_test_target} TYPE) + if (_sbml_test_type STREQUAL "EXECUTABLE") + target_precompile_headers(${_sbml_test_target} PRIVATE + + + + ) + endif() +endforeach() diff --git a/chaste_sbml/SbmlRefModels/test/TestSbmlMath.hpp b/chaste_sbml/SbmlRefModels/test/TestSbmlMath.hpp index f592cb63..daa40883 100644 --- a/chaste_sbml/SbmlRefModels/test/TestSbmlMath.hpp +++ b/chaste_sbml/SbmlRefModels/test/TestSbmlMath.hpp @@ -154,6 +154,14 @@ class TestSbmlMath : public CxxTest::TestSuite TS_ASSERT_EQUALS(sm::not_(true), false); TS_ASSERT_EQUALS(sm::not_(false), true); } + void TestImplies() + { + // implies(a, b) == (not a) or b + TS_ASSERT_EQUALS(sm::implies(true, true), true); + TS_ASSERT_EQUALS(sm::implies(true, false), false); + TS_ASSERT_EQUALS(sm::implies(false, true), true); + TS_ASSERT_EQUALS(sm::implies(false, false), true); + } void TestXor() { @@ -166,6 +174,11 @@ class TestSbmlMath : public CxxTest::TestSuite TS_ASSERT_EQUALS(sm::xor_(false, false), false); TS_ASSERT_EQUALS(sm::xor_(false, true, false), true); TS_ASSERT_EQUALS(sm::xor_(true, false, true), false); + + // Numeric-boolean (double) operands: cast to bool before xoring. + TS_ASSERT_EQUALS(sm::xor_(1.0, 0.0), true); + TS_ASSERT_EQUALS(sm::xor_(1.0, 1.0), false); + TS_ASSERT_EQUALS(sm::xor_(1.0, 0.0, 1.0), false); } // Relational ================================= @@ -361,6 +374,7 @@ class TestSbmlMath : public CxxTest::TestSuite void TestMax() { + TS_ASSERT_EQUALS(sm::max(2.0), 2.0); // single argument TS_ASSERT_EQUALS(sm::max(1.0, 2.0), 2.0); TS_ASSERT_EQUALS(sm::max(2.0, 1.0), 2.0); TS_ASSERT_EQUALS(sm::max(1.0, 2.0, 3.0), 3.0); @@ -374,6 +388,7 @@ class TestSbmlMath : public CxxTest::TestSuite void TestMin() { + TS_ASSERT_EQUALS(sm::min(2.0), 2.0); // single argument TS_ASSERT_EQUALS(sm::min(1.0, 2.0), 1.0); TS_ASSERT_EQUALS(sm::min(2.0, 1.0), 1.0); TS_ASSERT_EQUALS(sm::min(3.0, 2.0, 1.0), 1.0); @@ -394,6 +409,12 @@ class TestSbmlMath : public CxxTest::TestSuite TS_ASSERT_EQUALS(sm::piecewise(1, true, 2), 1); TS_ASSERT_EQUALS(sm::piecewise(1, false, 2), 2); + + // No otherwise clause: the value when the condition holds, NaN when it does not. + TS_ASSERT_EQUALS(sm::piecewise(1.0, true), 1.0); + TS_ASSERT(std::isnan(sm::piecewise(1.0, false))); + TS_ASSERT_EQUALS(sm::piecewise(1.0, true, 2.0, false), 1.0); + TS_ASSERT(std::isnan(sm::piecewise(1.0, false, 2.0, false))); } void TestQuotient() diff --git a/chaste_sbml/SbmlRefModels/test/cases/.gitignore b/chaste_sbml/SbmlRefModels/test/cases/.gitignore new file mode 100644 index 00000000..855b75d4 --- /dev/null +++ b/chaste_sbml/SbmlRefModels/test/cases/.gitignore @@ -0,0 +1,3 @@ +semantic/ +stochastic/ +syntactic/ diff --git a/chaste_sbml/SbmlRefModels/test/cases/README.md b/chaste_sbml/SbmlRefModels/test/cases/README.md new file mode 100644 index 00000000..1731ab44 --- /dev/null +++ b/chaste_sbml/SbmlRefModels/test/cases/README.md @@ -0,0 +1,10 @@ +# SBML Test Suite Cases + +## Semantic +Case support tracked in test/data/sbml_test_suite_status.csv + +## Syntactic +Covered by libSBML + +## Stochastic +Not supported \ No newline at end of file diff --git a/chaste_sbml/SbmlRefModels/test/data/sbml_test_suite_status.csv b/chaste_sbml/SbmlRefModels/test/data/sbml_test_suite_status.csv index 1447d775..ab3c1b76 100644 --- a/chaste_sbml/SbmlRefModels/test/data/sbml_test_suite_status.csv +++ b/chaste_sbml/SbmlRefModels/test/data/sbml_test_suite_status.csv @@ -7,11 +7,11 @@ test,status,notes 00006,pass,ok 00007,pass,ok 00008,pass,ok -00009,unsupported,no odes +00009,pass,ok 00010,pass,ok 00011,pass,ok 00012,pass,ok -00013,unsupported,no odes +00013,pass,ok 00014,pass,ok 00015,pass,ok 00016,pass,ok @@ -27,8 +27,8 @@ test,status,notes 00026,pass,ok 00027,pass,ok 00028,pass,ok -00029,unsupported,no odes -00030,unsupported,no odes +00029,pass,ok +00030,pass,ok 00031,pass,ok 00032,pass,ok 00033,pass,ok @@ -116,7 +116,7 @@ test,status,notes 00115,pass,ok 00116,pass,ok 00117,pass,ok -00118,unsupported,no odes +00118,pass,ok 00119,pass,ok 00120,pass,ok 00121,pass,ok @@ -172,7 +172,7 @@ test,status,notes 00171,pass,ok 00172,pass,ok 00173,pass,ok -00174,unsupported,no odes +00174,pass,ok 00175,pass,ok 00176,pass,ok 00177,pass,ok @@ -211,7 +211,7 @@ test,status,notes 00210,pass,ok 00211,pass,ok 00212,pass,ok -00213,unsupported,no odes +00213,pass,ok 00214,pass,ok 00215,pass,ok 00216,pass,ok @@ -226,7 +226,7 @@ test,status,notes 00225,pass,ok 00226,pass,ok 00227,pass,ok -00228,unsupported,no odes +00228,pass,ok 00229,pass,ok 00230,pass,ok 00231,pass,ok @@ -577,7 +577,7 @@ test,status,notes 00576,unsupported,algebraic rule 00577,pass,ok 00578,pass,ok -00579,unsupported,no odes +00579,pass,ok 00580,pass,ok 00581,pass,ok 00582,pass,ok @@ -918,54 +918,54 @@ test,status,notes 00917,pass,ok 00918,pass,ok 00919,pass,ok -00920,unsupported,no odes -00921,unsupported,no odes -00922,unsupported,no odes -00923,unsupported,no odes -00924,unsupported,no odes -00925,unsupported,no odes +00920,pass,ok +00921,pass,ok +00922,pass,ok +00923,pass,ok +00924,pass,ok +00925,pass,ok 00926,pass,ok 00927,pass,ok 00928,pass,ok 00929,pass,ok -00930,unsupported,no odes -00931,unsupported,no odes +00930,pass,ok +00931,pass,ok 00932,unsupported,delay 00933,unsupported,delay -00934,unsupported,no odes -00935,unsupported,no odes -00936,unsupported,no odes -00937,unsupported,no odes +00934,unsupported,event execution semantics +00935,unsupported,event execution semantics +00936,unsupported,delay +00937,unsupported,delay 00938,unsupported,delay 00939,unsupported,delay -00940,unsupported,no odes -00941,unsupported,no odes -00942,unsupported,no odes -00943,unsupported,no odes +00940,unsupported,delay +00941,unsupported,delay +00942,unsupported,delay +00943,unsupported,delay 00944,pass,ok 00945,pass,ok 00946,pass,ok 00947,pass,ok 00948,pass,ok -00949,unsupported,no odes -00950,unsupported,no odes -00951,unsupported,no odes -00952,unsupported,no odes -00953,unsupported,no odes -00954,unsupported,no odes -00955,unsupported,no odes -00956,unsupported,no odes -00957,unsupported,no odes -00958,unsupported,no odes -00959,unsupported,no odes -00960,unsupported,no odes -00961,unsupported,no odes -00962,unsupported,no odes -00963,unsupported,no odes -00964,unsupported,no odes -00965,unsupported,no odes -00966,unsupported,no odes -00967,unsupported,no odes +00949,pass,ok +00950,pass,ok +00951,pass,ok +00952,unsupported,random event execution +00953,unsupported,event execution semantics +00954,pass,ok +00955,pass,ok +00956,pass,ok +00957,pass,ok +00958,pass,ok +00959,pass,ok +00960,pass,ok +00961,pass,ok +00962,unsupported,random event execution +00963,unsupported,event execution semantics +00964,unsupported,random event execution +00965,unsupported,random event execution +00966,unsupported,random event execution +00967,unsupported,event execution semantics 00968,pass,ok 00969,pass,ok 00970,pass,ok @@ -976,9 +976,9 @@ test,status,notes 00975,pass,ok 00976,pass,ok 00977,pass,ok -00978,unsupported,no odes -00979,unsupported,no odes -00980,unsupported,no odes +00978,unsupported,event execution semantics +00979,pass,ok +00980,unsupported,delay 00981,unsupported,delay 00982,unsupported,delay 00983,unsupported,algebraic rule @@ -993,9 +993,9 @@ test,status,notes 00992,pass,ok 00993,unsupported,algebraic rule 00994,pass,ok -00995,unsupported,no odes +00995,pass,ok 00996,pass,ok -00997,unsupported,no odes +00997,pass,ok 00998,pass,ok 00999,pass,ok 01000,unsupported,delay @@ -1110,33 +1110,33 @@ test,status,notes 01109,pass,ok 01110,pass,ok 01111,pass,ok -01112,unsupported,no odes -01113,unsupported,no odes -01114,unsupported,no odes -01115,unsupported,no odes -01116,unsupported,no odes +01112,pass,ok +01113,pass,ok +01114,pass,ok +01115,pass,ok +01116,pass,ok 01117,pass,ok 01118,pass,ok -01119,unsupported,no odes +01119,unsupported,delay 01120,unsupported,delay 01121,pass,ok 01122,pass,ok 01123,pass,ok -01124,unsupported,no odes -01125,unsupported,no odes -01126,unsupported,no odes +01124,pass,ok +01125,pass,ok +01126,pass,ok 01127,pass,ok -01128,unsupported,no odes -01129,unsupported,no odes -01130,unsupported,no odes -01131,unsupported,no odes -01132,unsupported,no odes -01133,unsupported,no odes -01134,unsupported,no odes -01135,unsupported,no odes -01136,unsupported,no odes -01137,unsupported,no odes -01138,unsupported,no odes +01128,pass,ok +01129,pass,ok +01130,pass,ok +01131,pass,ok +01132,pass,ok +01133,pass,ok +01134,pass,ok +01135,pass,ok +01136,pass,ok +01137,pass,ok +01138,pass,ok 01139,pass,ok 01140,pass,ok 01141,pass,ok @@ -1144,23 +1144,23 @@ test,status,notes 01143,pass,ok 01144,pass,ok 01145,pass,ok -01146,unsupported,no odes -01147,unsupported,no odes +01146,pass,ok +01147,pass,ok 01148,pass,ok -01149,unsupported,no odes +01149,pass,ok 01150,pass,ok 01151,pass,ok -01152,unsupported,no odes -01153,unsupported,no odes -01154,unsupported,no odes -01155,unsupported,no odes -01156,unsupported,no odes -01157,unsupported,no odes -01158,unsupported,no odes +01152,pass,ok +01153,pass,ok +01154,pass,ok +01155,pass,ok +01156,unsupported,delay +01157,pass,ok +01158,pass,ok 01159,pass,ok 01160,pass,ok 01161,pass,ok -01162,unsupported,no odes +01162,pass,ok 01163,pass,ok 01164,pass,ok 01165,pass,ok @@ -1169,13 +1169,13 @@ test,status,notes 01168,pass,ok 01169,pass,ok 01170,pass,ok -01171,unsupported,no odes +01171,pass,ok 01172,pass,ok -01173,unsupported,no odes +01173,unsupported,delay 01174,unsupported,algebraic rule 01175,pass,ok 01176,unsupported,delay -01177,unsupported,no odes +01177,unsupported,delay 01178,pass,ok 01179,pass,ok 01180,pass,ok @@ -1195,7 +1195,7 @@ test,status,notes 01194,unsupported,flux balance 01195,unsupported,flux balance 01196,unsupported,flux balance -01197,unsupported,flux balance +01197,pass,ok 01198,pass,ok 01199,pass,ok 01200,pass,ok @@ -1207,16 +1207,16 @@ test,status,notes 01206,pass,ok 01207,pass,ok 01208,pass,ok -01209,unsupported,no odes -01210,unsupported,no odes -01211,unsupported,no odes -01212,unsupported,no odes -01213,unsupported,no odes -01214,unsupported,no odes +01209,pass,ok +01210,pass,ok +01211,pass,ok +01212,pass,ok +01213,unsupported,delay +01214,pass,ok 01215,pass,ok -01216,unsupported,no odes -01217,unsupported,no odes -01218,unsupported,no odes +01216,pass,ok +01217,pass,ok +01218,pass,ok 01219,pass,ok 01220,pass,ok 01221,pass,ok @@ -1232,22 +1232,22 @@ test,status,notes 01231,pass,ok 01232,pass,ok 01233,pass,ok -01234,unsupported,no odes -01235,unsupported,no odes -01236,unsupported,no odes -01237,unsupported,no odes -01238,unsupported,no odes -01239,unsupported,no odes -01240,unsupported,no odes -01241,unsupported,no odes -01242,unsupported,no odes -01243,unsupported,no odes +01234,pass,ok +01235,pass,ok +01236,pass,ok +01237,pass,ok +01238,pass,ok +01239,pass,ok +01240,pass,ok +01241,pass,ok +01242,pass,ok +01243,pass,ok 01244,unsupported,algebraic rule -01245,unsupported,no odes +01245,pass,ok 01246,pass,ok -01247,unsupported,no odes -01248,unsupported,no odes -01249,unsupported,no odes +01247,pass,ok +01248,pass,ok +01249,pass,ok 01250,pass,ok 01251,pass,ok 01252,pass,ok @@ -1269,27 +1269,27 @@ test,status,notes 01268,unsupported,delay 01269,unsupported,delay 01270,unsupported,delay -01271,unsupported,no odes -01272,unsupported,no odes -01273,unsupported,no odes -01274,unsupported,no odes -01275,unsupported,no odes -01276,unsupported,no odes -01277,unsupported,no odes -01278,unsupported,no odes -01279,unsupported,no odes -01280,unsupported,no odes -01281,unsupported,no odes -01282,unsupported,no odes -01283,unsupported,no odes -01284,unsupported,no odes -01285,unsupported,no odes -01286,unsupported,no odes -01287,unsupported,no odes +01271,pass,ok +01272,pass,ok +01273,pass,ok +01274,pass,ok +01275,pass,ok +01276,pass,ok +01277,pass,ok +01278,pass,ok +01279,pass,ok +01280,pass,ok +01281,pass,ok +01282,pass,ok +01283,pass,ok +01284,pass,ok +01285,pass,ok +01286,pass,ok +01287,unsupported,delay 01288,pass,ok -01289,unsupported,no odes +01289,pass,ok 01290,pass,ok -01291,unsupported,no odes +01291,pass,ok 01292,unsupported,algebraic rule 01293,pass,ok 01294,pass,ok @@ -1298,50 +1298,50 @@ test,status,notes 01297,pass,ok 01298,pass,ok 01299,unsupported,delay -01300,unsupported,no odes -01301,unsupported,no odes +01300,pass,ok +01301,pass,ok 01302,pass,ok -01303,unsupported,no odes -01304,unsupported,no odes -01305,unsupported,no odes -01306,unsupported,no odes +01303,pass,ok +01304,pass,ok +01305,unsupported,delay +01306,pass,ok 01307,pass,ok 01308,pass,ok 01309,pass,ok -01310,unsupported,no odes -01311,unsupported,no odes -01312,unsupported,no odes -01313,unsupported,no odes -01314,unsupported,no odes -01315,unsupported,no odes -01316,unsupported,no odes -01317,unsupported,no odes -01318,unsupported,no odes -01319,unsupported,no odes +01310,pass,ok +01311,pass,ok +01312,pass,ok +01313,pass,ok +01314,pass,ok +01315,pass,ok +01316,pass,ok +01317,pass,ok +01318,unsupported,delay +01319,unsupported,delay 01320,unsupported,delay 01321,pass,ok 01322,pass,ok -01323,unsupported,no odes +01323,pass,ok 01324,unsupported,delay 01325,unsupported,delay 01326,unsupported,delay 01327,unsupported,delay -01328,unsupported,no odes -01329,unsupported,no odes -01330,unsupported,no odes -01331,unsupported,no odes -01332,unsupported,no odes -01333,unsupported,no odes -01334,unsupported,no odes -01335,unsupported,no odes -01336,unsupported,no odes -01337,unsupported,no odes +01328,unsupported,delay +01329,unsupported,delay +01330,pass,ok +01331,unsupported,event execution semantics +01332,pass,ok +01333,unsupported,event execution semantics +01334,pass,ok +01335,unsupported,delay +01336,pass,ok +01337,unsupported,event execution semantics 01338,pass,ok 01339,pass,ok 01340,pass,ok -01341,unsupported,no odes -01342,unsupported,no odes -01343,unsupported,no odes +01341,pass,ok +01342,pass,ok +01343,pass,ok 01344,pass,ok 01345,pass,ok 01346,pass,ok @@ -1351,30 +1351,30 @@ test,status,notes 01350,unsupported,algebraic rule 01351,pass,ok 01352,pass,ok -01353,unsupported,no odes -01354,unsupported,no odes -01355,unsupported,no odes -01356,unsupported,no odes -01357,unsupported,no odes -01358,unsupported,no odes +01353,pass,ok +01354,pass,ok +01355,pass,ok +01356,pass,ok +01357,pass,ok +01358,unsupported,delay 01359,unsupported,algebraic rule 01360,pass,ok 01361,pass,ok -01362,unsupported,no odes -01363,unsupported,no odes -01364,unsupported,no odes -01365,unsupported,no odes -01366,unsupported,no odes -01367,unsupported,no odes +01362,pass,ok +01363,pass,ok +01364,pass,ok +01365,pass,ok +01366,pass,ok +01367,unsupported,delay 01368,unsupported,algebraic rule 01369,pass,ok 01370,pass,ok -01371,unsupported,no odes -01372,unsupported,no odes -01373,unsupported,no odes -01374,unsupported,no odes -01375,unsupported,no odes -01376,unsupported,no odes +01371,pass,ok +01372,pass,ok +01373,pass,ok +01374,pass,ok +01375,pass,ok +01376,unsupported,delay 01377,unsupported,algebraic rule 01378,pass,ok 01379,pass,ok @@ -1388,11 +1388,11 @@ test,status,notes 01387,pass,ok 01388,pass,ok 01389,pass,ok -01390,unsupported,no odes -01391,unsupported,no odes -01392,unsupported,no odes -01393,unsupported,no odes -01394,unsupported,no odes +01390,pass,ok +01391,pass,ok +01392,pass,ok +01393,pass,ok +01394,pass,ok 01395,pass,ok 01396,unsupported,fast reaction 01397,unsupported,fast reaction @@ -1414,7 +1414,7 @@ test,status,notes 01413,unsupported,delay 01414,unsupported,delay 01415,unsupported,delay -01416,unsupported,no odes +01416,unsupported,delay 01417,unsupported,delay 01418,unsupported,delay 01419,unsupported,delay @@ -1459,13 +1459,13 @@ test,status,notes 01458,pass,ok 01459,pass,ok 01460,pass,ok -01461,unsupported,no odes +01461,pass,ok 01462,pass,ok 01463,pass,ok 01464,pass,ok 01465,pass,ok -01466,unsupported,no odes -01467,unsupported,no odes +01466,pass,ok +01467,pass,ok 01468,pass,ok 01469,pass,ok 01470,pass,ok @@ -1483,19 +1483,19 @@ test,status,notes 01482,unsupported,algebraic rule 01483,unsupported,algebraic rule 01484,unsupported,algebraic rule -01485,unsupported,no odes -01486,unsupported,no odes -01487,unsupported,no odes -01488,unsupported,no odes -01489,unsupported,no odes -01490,unsupported,no odes -01491,unsupported,no odes +01485,pass,ok +01486,pass,ok +01487,pass,ok +01488,pass,ok +01489,pass,ok +01490,pass,ok +01491,pass,ok 01492,pass,ok 01493,pass,ok -01494,unsupported,no odes -01495,unsupported,no odes -01496,unsupported,no odes -01497,unsupported,no odes +01494,pass,ok +01495,pass,ok +01496,pass,ok +01497,pass,ok 01498,pass,ok 01499,unsupported,algebraic rule 01500,unsupported,algebraic rule @@ -1520,7 +1520,7 @@ test,status,notes 01519,unsupported,delay 01520,unsupported,delay 01521,unsupported,delay -01522,unsupported,no odes +01522,unsupported,delay 01523,unsupported,delay 01524,unsupported,delay 01525,unsupported,delay @@ -1528,10 +1528,10 @@ test,status,notes 01527,pass,ok 01528,unsupported,delay 01529,unsupported,delay -01530,unsupported,no odes -01531,partly_supported,clustered events +01530,pass,ok +01531,unsupported,event execution semantics 01532,unsupported,delay -01533,unsupported,no odes +01533,pass,ok 01534,unsupported,delay 01535,unsupported,delay 01536,unsupported,delay @@ -1586,23 +1586,23 @@ test,status,notes 01585,unsupported,delay 01586,unsupported,delay 01587,unsupported,delay -01588,unsupported,no odes +01588,unsupported,random event execution 01589,unsupported,algebraic rule -01590,unsupported,no odes +01590,unsupported,delay 01591,unsupported,random event execution 01592,unsupported,delay 01593,unsupported,delay -01594,unsupported,no odes -01595,unsupported,no odes -01596,unsupported,no odes -01597,unsupported,no odes -01598,unsupported,no odes +01594,unsupported,delay +01595,unsupported,delay +01596,pass,ok +01597,unsupported,delay +01598,unsupported,delay 01599,unsupported,random event execution -01600,unsupported,no odes -01601,unsupported,no odes -01602,unsupported,no odes -01603,unsupported,no odes -01604,unsupported,no odes +01600,unsupported,delay +01601,unsupported,delay +01602,unsupported,delay +01603,unsupported,delay +01604,unsupported,delay 01605,unsupported,random event execution 01606,unsupported,flux balance 01607,unsupported,flux balance @@ -1633,13 +1633,13 @@ test,status,notes 01632,pass,ok 01633,pass,ok 01634,pass,ok -01635,unsupported,no odes +01635,pass,ok 01636,pass,ok 01637,pass,ok -01638,unsupported,no odes +01638,pass,ok 01639,pass,ok 01640,pass,ok -01641,unsupported,no odes +01641,pass,ok 01642,pass,ok 01643,pass,ok 01644,pass,ok @@ -1656,13 +1656,13 @@ test,status,notes 01655,pass,ok 01656,pass,ok 01657,pass,ok -01658,unsupported,no odes -01659,unsupported,no odes -01660,unsupported,no odes -01661,unsupported,no odes -01662,unsupported,no odes -01663,unsupported,no odes -01664,unsupported,no odes +01658,pass,ok +01659,unsupported,delay +01660,unsupported,delay +01661,unsupported,delay +01662,pass,ok +01663,pass,ok +01664,pass,ok 01665,pass,ok 01666,pass,ok 01667,pass,ok @@ -1691,14 +1691,14 @@ test,status,notes 01690,unsupported,delay 01691,unsupported,delay 01692,unsupported,delay -01693,unsupported,no odes -01694,unsupported,no odes -01695,unsupported,no odes -01696,unsupported,no odes -01697,unsupported,no odes -01698,unsupported,no odes -01699,unsupported,no odes -01700,unsupported,no odes +01693,pass,ok +01694,pass,ok +01695,pass,ok +01696,pass,ok +01697,pass,ok +01698,pass,ok +01699,pass,ok +01700,pass,ok 01701,unsupported,delay 01702,unsupported,delay 01703,unsupported,delay @@ -1752,12 +1752,12 @@ test,status,notes 01751,pass,ok 01752,pass,ok 01753,pass,ok -01754,unsupported,no odes -01755,unsupported,no odes -01756,unsupported,no odes -01757,unsupported,no odes -01758,unsupported,no odes -01759,unsupported,no odes +01754,unsupported,delay +01755,unsupported,delay +01756,unsupported,delay +01757,unsupported,delay +01758,unsupported,delay +01759,unsupported,delay 01760,pass,ok 01761,pass,ok 01762,pass,ok @@ -1777,10 +1777,10 @@ test,status,notes 01776,pass,ok 01777,pass,ok 01778,pass,ok -01779,unsupported,no odes -01780,unsupported,no odes -01781,unsupported,no odes -01782,unsupported,no odes +01779,pass,ok +01780,pass,ok +01781,pass,ok +01782,pass,ok 01783,pass,ok 01784,pass,ok 01785,unsupported,algebraic rule @@ -1794,29 +1794,29 @@ test,status,notes 01793,unsupported,algebraic rule 01794,unsupported,algebraic rule 01795,pass,ok -01796,unsupported,no odes +01796,pass,ok 01797,pass,ok 01798,unsupported,delay 01799,pass,ok 01800,pass,ok 01801,pass,ok 01802,pass,ok -01803,unsupported,no odes -01804,unsupported,no odes -01805,unsupported,no odes -01806,unsupported,no odes +01803,pass,ok +01804,pass,ok +01805,pass,ok +01806,pass,ok 01807,pass,ok 01808,pass,ok 01809,pass,ok -01810,unsupported,no odes -01811,unsupported,no odes -01812,unsupported,no odes -01813,unsupported,no odes -01814,unsupported,no odes -01815,unsupported,no odes -01816,unsupported,no odes -01817,unsupported,no odes -01818,unsupported,no odes -01819,unsupported,no odes -01820,unsupported,no odes -01821,unsupported,no odes +01810,pass,ok +01811,pass,ok +01812,pass,ok +01813,pass,ok +01814,pass,ok +01815,pass,ok +01816,pass,ok +01817,pass,ok +01818,pass,ok +01819,pass,ok +01820,pass,ok +01821,pass,ok diff --git a/chaste_sbml/SbmlRefModels/test/reference/TestChen2004SbmlOdeSystem.hpp b/chaste_sbml/SbmlRefModels/test/reference/TestChen2004SbmlOdeSystem.hpp index 3e932785..3d62b88f 100644 --- a/chaste_sbml/SbmlRefModels/test/reference/TestChen2004SbmlOdeSystem.hpp +++ b/chaste_sbml/SbmlRefModels/test/reference/TestChen2004SbmlOdeSystem.hpp @@ -66,7 +66,7 @@ class TestChen2004SbmlOdeSystem : public AbstractCellBasedTestSuite { private: const unsigned ODE_SIZE = 36u; - const unsigned NUM_DERIVED_QUANTITIES = 187u; + const unsigned NUM_DERIVED_QUANTITIES = 184u; std::vector default_initial_conditions = { 0.008473, // BUD @@ -470,7 +470,6 @@ class TestChen2004SbmlOdeSystem : public AbstractCellBasedTestSuite std::vector dq_names = { "cell", "BCK2", - "BUB2", "CDC14T", "CDC15i", "CDC6T", @@ -479,8 +478,6 @@ class TestChen2004SbmlOdeSystem : public AbstractCellBasedTestSuite "CLB5T", "CLN3", "IE", - "LTE1", - "MAD2", "MCM1", "NET1T", "PE", @@ -524,7 +521,6 @@ class TestChen2004SbmlOdeSystem : public AbstractCellBasedTestSuite std::vector dqs_expected = { 1.0, // cell 0.065125026, // BCK2 - 0.2, // BUB2 2.117884, // CDC14T 0.34346699999999997, // CDC15i 0.3866693, // CDC6T @@ -533,8 +529,6 @@ class TestChen2004SbmlOdeSystem : public AbstractCellBasedTestSuite 0.1289119, // CLB5T 0.06694509131879892, // CLN3 0.8985, // IE - 0.1, // LTE1 - 0.01, // MAD2 0.4690076182110798, // MCM1 2.638456, // NET1T 0.6986870000000001, // PE diff --git a/chaste_sbml/_config.py b/chaste_sbml/_config.py index e9d70ebd..cc37b5fb 100644 --- a/chaste_sbml/_config.py +++ b/chaste_sbml/_config.py @@ -16,6 +16,10 @@ DERIVATIVE_PREFIX = "d_" DERIVATIVE_SUFFIX = "_dt" +# Base id for the synthetic zero-derivative state variable added to models with no ODEs, so the +# generated ODE system always has at least one variable for the solver to integrate. +PLACEHOLDER_STATE_ID = "placeholder" + class EquationType(Enum): """Enumeration of equation types for code generation.""" @@ -47,6 +51,19 @@ class VarType(Enum): UNKNOWN = 5 +class DerivedQuantityKind(Enum): + """Kind of derived quantity, which determines how it is declared and computed. + + NORMAL is declared as a member and computed from an equation; REACTION is a reaction flux + (declared and computed by the reaction machinery, only exposed here as an output); CONVERSION + is an amount/concentration conversion computed in ComputeDerivedQuantities. + """ + + NORMAL = 0 + REACTION = 1 + CONVERSION = 2 + + class ModelType(Enum): """Enumeration of model types for code generation.""" diff --git a/chaste_sbml/_expressions.py b/chaste_sbml/_expressions.py index e3a304fe..06faecad 100644 --- a/chaste_sbml/_expressions.py +++ b/chaste_sbml/_expressions.py @@ -6,10 +6,29 @@ into the model. """ +import math import re from typing import TYPE_CHECKING, Optional -from libsbml import AST_FUNCTION_DELAY, AST_NAME, AST_NAME_AVOGADRO, formulaToL3String +from libsbml import ( + AST_CONSTANT_E, + AST_CONSTANT_PI, + AST_FUNCTION, + AST_FUNCTION_DELAY, + AST_FUNCTION_POWER, + AST_FUNCTION_RATE_OF, + AST_NAME, + AST_NAME_AVOGADRO, + AST_NAME_TIME, + AST_POWER, + AST_REAL, + AST_RELATIONAL_EQ, + AST_RELATIONAL_GEQ, + AST_RELATIONAL_GT, + AST_RELATIONAL_LEQ, + AST_RELATIONAL_LT, + formulaToL3String, +) from ._config import CHASTE_PREFIX, PREFIX_SEP, VarType @@ -55,21 +74,89 @@ def substitute_ast_names(node: "ASTNode", replacements: dict) -> "ASTNode": return result -def replace_avogadro_csymbol(node: "ASTNode", placeholder: str) -> None: - """Recursively rename avogadro csymbol nodes to a placeholder identifier. - - The avogadro csymbol and a parameter both named 'avogadro' are distinct AST nodes that - ``formulaToL3String`` renders identically. Renaming the csymbol to a placeholder keeps it - distinct so it can be mapped to ``sm::AVOGADRO`` while the parameter keeps its own name. +# Placeholder ids the SBML math symbols are renamed to. A parameter may share a symbol's spelling +# (e.g. a parameter named 'pi' or 'time'), and ``formulaToL3String`` renders the symbol and the +# parameter identically; renaming the symbol node to a unique placeholder keeps them distinct, so +# the symbol maps to its C++ value (below) while the parameter keeps its own name. +AVOGADRO_PLACEHOLDER = f"{CHASTE_PREFIX}{PREFIX_SEP}avogadro" +PI_PLACEHOLDER = f"{CHASTE_PREFIX}{PREFIX_SEP}pi" +E_PLACEHOLDER = f"{CHASTE_PREFIX}{PREFIX_SEP}exponentiale" +TIME_PLACEHOLDER = f"{CHASTE_PREFIX}{PREFIX_SEP}time" +INFINITY_PLACEHOLDER = f"{CHASTE_PREFIX}{PREFIX_SEP}infinity" +NEG_INFINITY_PLACEHOLDER = f"{CHASTE_PREFIX}{PREFIX_SEP}neg_infinity" +NAN_PLACEHOLDER = f"{CHASTE_PREFIX}{PREFIX_SEP}nan" + +_SYMBOL_PLACEHOLDERS = { + AST_NAME_AVOGADRO: AVOGADRO_PLACEHOLDER, + AST_CONSTANT_PI: PI_PLACEHOLDER, + AST_CONSTANT_E: E_PLACEHOLDER, + AST_NAME_TIME: TIME_PLACEHOLDER, +} + + +def replace_constant_symbols(node: "ASTNode") -> None: + """Recursively rename SBML math-symbol nodes (avogadro, pi, exponentiale, time, inf, NaN). + + Each of these symbols renders identically to a same-named parameter, so it is renamed to a + unique placeholder that maps to its C++ value, keeping it distinct from a parameter of the same + spelling. Handling this at the AST level (rather than string-matching the rendered name) means a + real variable named ``t``, ``pi`` etc. is never mistaken for the symbol. The MathML constants + ```` and ```` are especially treacherous: libsbml stores them as ordinary + real nodes that render as ``INF``/``NaN``, indistinguishable from a reference to a parameter so + named (see test-suite cases 1811/1813), so they are matched here by their non-finite value. :param node: The root AST node. - :param placeholder: The identifier to rename avogadro csymbol nodes to. """ - if node.getType() == AST_NAME_AVOGADRO: + placeholder = _SYMBOL_PLACEHOLDERS.get(node.getType()) + if placeholder is not None: node.setType(AST_NAME) node.setName(placeholder) + elif node.getType() == AST_REAL: + value = node.getReal() + if math.isnan(value): + node.setType(AST_NAME) + node.setName(NAN_PLACEHOLDER) + elif math.isinf(value): + node.setType(AST_NAME) + node.setName(NEG_INFINITY_PLACEHOLDER if value < 0 else INFINITY_PLACEHOLDER) for i in range(node.getNumChildren()): - replace_avogadro_csymbol(node.getChild(i), placeholder) + replace_constant_symbols(node.getChild(i)) + + +def resolve_rate_of( + node: "ASTNode", + variable_types: dict, + state_variables: list["StateVariable"], + local_param_ids: set, +) -> None: + """Recursively resolve ``rateOf(x)`` nodes to a state variable's derivative id, or ``0``. + + ``rateOf`` of a state variable is that variable's time derivative (its ``d__dt`` id); of + anything constant -- a parameter, or a kinetic-law local parameter that shadows a global -- it + is zero. Resolving on the AST (rather than by regex on the rendered string) is robust to the + argument's spelling and to surrounding parentheses. + + :param node: The root AST node. + :param variable_types: Mapping of model id to :class:`VarType`. + :param state_variables: The model's state variables (to look up a derivative id). + :param local_param_ids: Ids of local parameters in scope (constant, so their rate is zero). + """ + if node.getType() == AST_FUNCTION_RATE_OF and node.getNumChildren() >= 1: + var = node.getChild(0).getName() or "" + derivative = None + if var not in local_param_ids and variable_types.get(var, VarType.UNKNOWN) == VarType.STATE_VARIABLE: + derivative = next((sv.derivative_id for sv in state_variables if sv.id == var), None) + while node.getNumChildren() > 0: + node.removeChild(0) + if derivative is not None: + node.setType(AST_NAME) + node.setName(derivative) + else: + node.setType(AST_REAL) + node.setValue(0.0) + return + for i in range(node.getNumChildren()): + resolve_rate_of(node.getChild(i), variable_types, state_variables, local_param_ids) def strip_ast_units(node: "ASTNode") -> None: @@ -77,16 +164,52 @@ def strip_ast_units(node: "ASTNode") -> None: SBML allows numeric literals to carry units annotations (e.g. ````). ``formulaToL3String`` includes these in its output (e.g. ``0.00015 mole``), which is not - valid C++. Stripping them here is safe because the units carry no mathematical information. + valid C++. Stripping them here is safe because the units are not used. :param node: The root AST node. """ + # TODO: Handle units properly if node.isSetUnits(): node.unsetUnits() for i in range(node.getNumChildren()): strip_ast_units(node.getChild(i)) +# SBML MathML relational operators, and the SbmlMath (``sm::``) helper each maps to. +_NARY_RELATIONAL_NAMES = { + AST_RELATIONAL_LT: "lt", + AST_RELATIONAL_GT: "gt", + AST_RELATIONAL_LEQ: "leq", + AST_RELATIONAL_GEQ: "geq", + AST_RELATIONAL_EQ: "eq", +} + + +def rewrite_nary_relational(node: Optional["ASTNode"]) -> None: + """Rewrite n-ary (>2 operand) relational AST nodes into function-call form, in place. + + SBML MathML relationals are n-ary: ``lt(a, b, c)`` means ``a < b < c`` + (i.e. ``a < b`` and ``b < c``), and ``eq(a, b, c)`` means all equal. However, + ``formulaToL3String`` renders them as a chained infix expression (``a < b < c``), + which C++ mis-evaluates left-to-right as ``(a < b) < c``. Converting each + such node to a function call (rendered ``lt(a, b, c)``, later mapped to the + variadic ``sm::lt`` helper that implements the correct semantics) fixes this. + Binary relationals are left as infix operators, so existing generated code + is unchanged. ``neq`` is binary-only in MathML, so it is never n-ary and is + not rewritten here. + + :param node: The root AST node to rewrite. + """ + if node is None: + return + name = _NARY_RELATIONAL_NAMES.get(node.getType()) + if name is not None and node.getNumChildren() > 2: + node.setType(AST_FUNCTION) + node.setName(name) + for i in range(node.getNumChildren()): + rewrite_nary_relational(node.getChild(i)) + + def search_ast_type(root: Optional["ASTNode"], node_type: int) -> bool: """Recursively search the AST for a node of a certain type. @@ -109,126 +232,124 @@ def search_ast_type(root: Optional["ASTNode"], node_type: int) -> bool: return False -def convert_infix_operator_to_function_syntax(formula: str, operator: str, function_name: str) -> str: - """Convert infix operator expressions to function syntax. +def rewrite_power(node: Optional["ASTNode"]) -> None: + """Rewrite power-operator (``a ^ b``) AST nodes into ``pow()`` function-call form, in place. - Example: with operator='^' and function_name='pow', rewrites - ``a ^ b`` to ``pow(a, b)``. + ``formulaToL3String`` renders a power as the infix ``a ^ b``, which is not + valid C++ (``^`` is bitwise xor). Converting each power node to a ``pow`` + call before stringifying (rendered ``pow(a, b)``, later mapped to ``std::pow`` + by the function-name token map) handles nesting correctly. Both power node + types are handled: ``AST_POWER`` (from ``^``/``pow`` in L3 infix, and from + MathML ```` that libsbml normalises) and ``AST_FUNCTION_POWER`` + (from other MathML ```` forms). - This parser handles parenthesized operands (including nested parentheses) - and simple symbolic/numeric tokens. - - :param formula: Formula text in SBML infix style. - :param operator: Infix operator token to rewrite. - :param function_name: Function name used for the replacement. - :return: Formula text with infix operators rewritten as function calls. + :param node: The root AST node to rewrite. """ - if not operator: - raise ValueError("operator must be a non-empty string") - if not function_name: - raise ValueError("function_name must be a non-empty string") - - token_chars = set("ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789_:.eE") - - def read_left_operand(s: str, operator_index: int) -> tuple[int, int] | None: - i = operator_index - 1 - while i >= 0 and s[i].isspace(): - i -= 1 - if i < 0: - return None - - if s[i] == ")": - depth = 1 - j = i - 1 - while j >= 0: - if s[j] == ")": - depth += 1 - elif s[j] == "(": - depth -= 1 - if depth == 0: - break - j -= 1 - if depth != 0: - return None - # j points to the matching '('; also include any function name before it - k = j - 1 - while k >= 0 and s[k] in token_chars: - k -= 1 - return (k + 1, i + 1) - - j = i - while j >= 0 and s[j] in token_chars: - j -= 1 - - start = j + 1 - end = i + 1 - if start >= end: - return None - return (start, end) - - def read_right_operand(s: str, operator_index: int) -> tuple[int, int] | None: - i = operator_index + len(operator) - n = len(s) - while i < n and s[i].isspace(): - i += 1 - if i >= n: - return None - - if s[i] == "(": - depth = 1 - j = i + 1 - while j < n: - if s[j] == "(": - depth += 1 - elif s[j] == ")": - depth -= 1 - if depth == 0: - return (i, j + 1) - j += 1 - return None - - j = i - if s[j] in "+-": - j += 1 - - while j < n and s[j] in token_chars: - j += 1 - - if j <= i: - return None - return (i, j) - - max_rewrites = 500 - rewrites = 0 - search_start = 0 - while rewrites < max_rewrites: - operator_index = formula.find(operator, search_start) - if operator_index < 0: - break - - left = read_left_operand(formula, operator_index) - right = read_right_operand(formula, operator_index) - if left is None or right is None: - search_start = operator_index + len(operator) - continue - - l_start, l_end = left - r_start, r_end = right - - lhs = formula[l_start:l_end].strip() - rhs = formula[r_start:r_end].strip() - replacement = f"{function_name}({lhs}, {rhs})" - - formula = formula[:l_start] + replacement + formula[r_end:] - search_start = l_start + len(replacement) - rewrites += 1 - - return formula - - -def convert_infix_power_to_function_syntax(formula: str) -> str: - """Convert infix power expressions (a ^ b) to function syntax pow(a, b).""" - return convert_infix_operator_to_function_syntax(formula=formula, operator="^", function_name="pow") + if node is None: + return + if node.getType() in (AST_POWER, AST_FUNCTION_POWER): + node.setType(AST_FUNCTION) + node.setName("pow") + for i in range(node.getNumChildren()): + rewrite_power(node.getChild(i)) + + +# SBML symbolic constants replaced with their C++ equivalents, keyed by the placeholder each symbol +# is renamed to (see replace_constant_symbols). Every one is matched at the AST level and renamed to +# its placeholder, so a same-named parameter (e.g. one literally called pi, inf or NaN) is never +# mistaken for the constant. (true/false deliberately absent.) +_CONSTANTS = { + AVOGADRO_PLACEHOLDER: "sm::AVOGADRO", + PI_PLACEHOLDER: "M_PI", + E_PLACEHOLDER: "M_E", + TIME_PLACEHOLDER: "time", + INFINITY_PLACEHOLDER: "std::numeric_limits::infinity()", + NEG_INFINITY_PLACEHOLDER: "-std::numeric_limits::infinity()", + NAN_PLACEHOLDER: "std::numeric_limits::quiet_NaN()", +} + +# SBML functions whose name matches the C++ (std::) equivalent. +_UNCHANGED_FUNCTIONS = { + "acos", + "acosh", + "asin", + "asinh", + "atan", + "atanh", + "ceil", + "cos", + "cosh", + "exp", + "floor", + "pow", + "sin", + "sinh", + "sqrt", + "tan", + "tanh", +} + +# SBML functions with a different name in C++ (std::). +_RENAMED_FUNCTIONS = { + "abs": "fabs", + "arccos": "acos", + "arccosh": "acosh", + "arcsin": "asin", + "arcsinh": "asinh", + "arctan": "atan", + "arctanh": "atanh", + "ceiling": "ceil", + "ln": "log", + "power": "pow", + "rem": "fmod", +} + +# SBML functions with custom implementations in SbmlMath (sm::). +_CUSTOM_FUNCTIONS = { + "and": "and_", + "acot": "acot", + "acoth": "acoth", + "acsc": "acsc", + "acsch": "acsch", + "asec": "asec", + "asech": "asech", + "arccot": "acot", + "arccoth": "acoth", + "arccsc": "acsc", + "arccsch": "acsch", + "arcsec": "asec", + "arcsech": "asech", + "cot": "cot", + "coth": "coth", + "csc": "csc", + "csch": "csch", + "eq": "eq", + "factorial": "factorial", + "geq": "geq", + "gt": "gt", + "implies": "implies", + "leq": "leq", + "log": "log", + "lt": "lt", + "max": "max", + "min": "min", + "neq": "neq", + "not": "not_", + "or": "or_", + "piecewise": "piecewise", + "plus": "plus", + "quotient": "quotient", + "root": "root", + "sec": "sec", + "sech": "sech", + "sqr": "sqr", + "times": "times", + "xor": "xor_", +} + +# TODO: From SBML Level 3 upwards, log defaults to base 10. +# SBML versions lower than 3 default to base e. +# See https://sbml.org/software/libsbml/5.18.0/docs/formatted/python-api/namespacelibsbml.html#a8e96a5a70569ae32655c6302638f6dc3 # noqa: B950 def formula_to_string( @@ -239,173 +360,62 @@ def formula_to_string( ) -> str: """Convert an AST math formula to an equivalent C++ string. - :param math: The AST math formula. - :param variable_types: Mapping of model id to :class:`VarType`, used to tell model variables - apart from SBML constants/functions of the same spelling. - :param state_variables: The model's state variables, used to resolve ``rateOf`` to a state - variable's derivative id. - :param local_parameters: Local parameters in scope (e.g. a reaction's kinetic-law - parameters). These shadow global symbols of the same name and are constant, so - ``rateOf`` applied to one is zero. + :param math: The AST math formula. It is not mutated. + :param variable_types: Mapping of model id to :class:`VarType`, used to tell + model variables apart from SBML constants/functions of the same spelling. + :param state_variables: The model's state variables, used to resolve ``rateOf`` + to a state variable's derivative id. + :param local_parameters: Local parameters in scope (e.g. a reaction's + kinetic-law parameters). These shadow global symbols of the same name + and are constant, so ``rateOf`` applied to one is zero. :return: The equivalent C++ string. """ - unsupported_functions = ["delay"] - for func in unsupported_functions: - if search_ast_type(math, AST_FUNCTION_DELAY): - raise NotImplementedError(f"SBML function not supported: '{func}'.") + if search_ast_type(math, AST_FUNCTION_DELAY): + raise NotImplementedError("SBML function not supported: 'delay'.") + + # Make a copy so the AST is not mutated. + math = math.deepCopy() + local_param_ids = {param.id for param in (local_parameters or [])} strip_ast_units(math) - # The avogadro csymbol and a parameter both named 'avogadro' are distinct AST nodes that - # render identically. Rename the csymbol to a placeholder (mapped to sm::AVOGADRO in the - # constants below) so it stays distinct from a same-named parameter. - avogadro_placeholder = f"{CHASTE_PREFIX}{PREFIX_SEP}avogadro" - replace_avogadro_csymbol(math, avogadro_placeholder) + # Modify the AST, renaming the avogadro/pi/exponentiale/time symbol nodes to + # placeholders (so each stays distinct from a same-named parameter), and + # resolve rateOf to a derivative id (or 0). + replace_constant_symbols(math) + resolve_rate_of(math, variable_types, state_variables, local_param_ids) + # Route operators that formulaToL3String would render as C++-invalid infix + # to function-call form: n-ary relationals (a < b < c) to the SBMLMath + # helpers (sm::lt(a, b, c)), and the power operator (a ^ b) to std::pow(). + # Binary relationals stay as infix operators (a < b). + rewrite_nary_relational(math) + rewrite_power(math) formula = formulaToL3String(math) - # Convert all integer literals to doubles to fix integer division. - # TODO: Perhaps instead of regex, traverse AST and convert some AST_INTEGER - # nodes to AST_REAL. This should only need to apply to numbers used in a division. - formula = re.sub(r"(?::infinity()", - "infinity": "std::numeric_limits::infinity()", - "nan": "NAN", - "notanumber": "NAN", - "pi": "M_PI", - "time": "time", - "t": "time", - "s": "time", - } - # skip: "true", "false" - - # SBML functions with same name as C++ equivalents - unchanged_functions = { - "acos", - "acosh", - "asin", - "asinh", - "atan", - "atanh", - "ceil", - "cos", - "cosh", - "exp", - "floor", - "pow", - "sin", - "sinh", - "sqrt", - "tan", - "tanh", - } - - # SBML functions with different names in C++ - renamed_functions = { - "abs": "fabs", - "arccos": "acos", - "arccosh": "acosh", - "arcsin": "asin", - "arcsinh": "asinh", - "arctan": "atan", - "arctanh": "atanh", - "ceiling": "ceil", - "ln": "log", - "power": "pow", - "rem": "fmod", - } - - # SBML functions with custom implementations - custom_functions = { - "and": "and_", - "acot": "acot", - "acoth": "acoth", - "acsc": "acsc", - "acsch": "acsch", - "asec": "asec", - "asech": "asech", - "arccot": "acot", - "arccoth": "acoth", - "arccsc": "acsc", - "arccsch": "acsch", - "arcsec": "asec", - "arcsech": "asech", - "cot": "cot", - "coth": "coth", - "csc": "csc", - "csch": "csch", - "eq": "eq", - "factorial": "factorial", - "geq": "geq", - "gt": "gt", - "leq": "leq", - "log": "log", - "lt": "lt", - "max": "max", - "min": "min", - "neq": "neq", - "not": "not_", - "or": "or_", - "piecewise": "piecewise", - "plus": "plus", - "quotient": "quotient", - "root": "root", - "sec": "sec", - "sech": "sech", - "sqr": "sqr", - "times": "times", - "xor": "xor_", - } - - # TODO: From SBML Level 3 upwards, log defaults to base 10. - # SBML versions lower than 3 default to base e. - # See https://sbml.org/software/libsbml/5.18.0/docs/formatted/python-api/namespacelibsbml.html#a8e96a5a70569ae32655c6302638f6dc3 # noqa: B950 + # Convert all integer literals to doubles to fix integer division. The lookbehinds + # keep the digits of a scientific-notation exponent (e.g. 1e-5, 1e+5) intact. + formula = re.sub(r"(? None: amt_rhs = f"{species_id} * {compartment_id}" amt_math = parseL3Formula(amt_rhs) - self._add_derived_quantity(amt_id, amt_label, None, amt_units, is_conversion=True) + self._add_derived_quantity(amt_id, amt_label, None, amt_units, kind=DerivedQuantityKind.CONVERSION) self._add_equation(var=amt_id, math=amt_math, eq_type=EquationType.CONVERSION) def _add_concentration(self, species: "Species") -> None: @@ -135,7 +137,7 @@ def _add_concentration(self, species: "Species") -> None: conc_rhs = f"{species_id} / {compartment_id}" conc_math = parseL3Formula(conc_rhs) - self._add_derived_quantity(conc_id, conc_label, None, conc_units, is_conversion=True) + self._add_derived_quantity(conc_id, conc_label, None, conc_units, kind=DerivedQuantityKind.CONVERSION) self._add_equation(var=conc_id, math=conc_math, eq_type=EquationType.CONVERSION) def _add_assignment_rule(self, id_: str, label: str, var: str, math: "ASTNode") -> None: @@ -154,8 +156,7 @@ def _add_derived_quantity( label: str, initial_value: Optional[float], units: str = NON_DIM_UNITS, - is_conversion: bool = False, - is_reaction: bool = False, + kind: DerivedQuantityKind = DerivedQuantityKind.NORMAL, ) -> None: """Add a derived quantity to the template variables. @@ -163,11 +164,10 @@ def _add_derived_quantity( :param label: The variable description. :param initial_value: The variable initial value. :param units: The variable units. - :param is_conversion: True if the derived quantity is an amount/concentration - conversion (computed in ComputeDerivedQuantities rather than stored as a member). - :param is_reaction: True if the derived quantity is a reaction flux exposed as an - output. Its member is already declared by the reaction loop, so it is not - re-declared, and its variable type stays VarType.REACTION. + :param kind: The kind of derived quantity. CONVERSION is an amount/concentration conversion + (computed in ComputeDerivedQuantities rather than stored as a member); REACTION is a + reaction flux exposed as an output -- its member is already declared by the reaction + loop, so it is not re-declared and its variable type stays VarType.REACTION. """ self._derived_quantities.append( DerivedQuantity( @@ -175,13 +175,11 @@ def _add_derived_quantity( label=label, index=len(self._derived_quantities), initial_value=initial_value, - is_conversion=is_conversion, - is_reaction=is_reaction, units=units, + kind=kind, + name=self._names.sbml_name(id_), ) ) - if not is_reaction: - self._variable_types[id_] = VarType.DERIVED_QUANTITY def _add_equation( self, @@ -243,7 +241,6 @@ def _add_function(self, id_: str, label: str, args: str, body: str) -> None: :param body: The function body. """ self._functions.append(Function(id=id_, label=label, index=len(self._functions), args=args, body=body)) - self._variable_types[id_] = VarType.FUNCTION def _add_initial_assignment(self, id_: str, label: str, var: str, math: Optional["ASTNode"] = None) -> None: """Add an initial assignment to the template variables. @@ -279,9 +276,9 @@ def _add_parameter( label=label, initial_value=initial_value, units=units, + name=self._names.sbml_name(id_), ) ) - self._variable_types[id_] = VarType.PARAMETER def _add_rate_rule(self, id_: str, label: str, var: str, math: "ASTNode") -> None: """Add a rate rule to the template variables. @@ -300,7 +297,6 @@ def _add_reaction(self, id_: str, label: str) -> None: :param label: The variable description. """ self._reactions.append(Reaction(index=len(self._reactions), id=id_, label=label)) - self._variable_types[id_] = VarType.REACTION def _add_state_variable( self, id_: str, label: str, initial_value: Optional[float], units: str = NON_DIM_UNITS @@ -320,13 +316,28 @@ def _add_state_variable( label=label, initial_value=initial_value, units=units, + name=self._names.sbml_name(id_), ) self._state_variables.append(state_var) - self._variable_types[id_] = VarType.STATE_VARIABLE return state_var + def _add_placeholder_state_variable(self) -> None: + """Add a synthetic state variable with a zero derivative for a model with no ODEs. + + Chaste's ODE system integrates a vector of state variables, so a model with no continuous + dynamics still needs at least one. The placeholder is constant (``dy/dt = 0``), starts at + zero and is not an SBML output, so it never appears in the expected results; it exists only + to give the solver a trivial state to advance in time while the model's real outputs are + recomputed each step. + """ + placeholder_id = self._names.reserve(PREFIX_SEP.join([CHASTE_PREFIX, PLACEHOLDER_STATE_ID])) + state_var = self._add_state_variable(placeholder_id, "Placeholder state variable (model has no ODEs)", 0.0) + # Initialise the member to zero + self._add_equation(var=placeholder_id, math=parseL3Formula("0.0"), eq_type=EquationType.INITIAL_VALUE) + self._add_equation(var=state_var.derivative_id, math=parseL3Formula("0.0"), eq_type=EquationType.DERIVATIVE) + def _add_stoichiometry_variable(self, species_reference: "SpeciesReference") -> None: """Add a stoichiometry variable to the template variables. @@ -449,9 +460,9 @@ def _update_ode(species_ref: "SpeciesReference", rxn: "SbmlReaction", add: bool) math = rule.math odes[var] = math - if not odes: - raise NotImplementedError("Models without ODEs are not supported.") - + # An empty odes dict is allowed: a model with no reactions or rate rules has no + # continuous dynamics, and build() synthesises a placeholder state variable so the + # generated ODE system still has something for the solver to integrate. self._odes = odes def _total_time_derivative(self, ast_node: "ASTNode") -> str: @@ -531,12 +542,16 @@ def _format_equations(self) -> None: def _format_events(self) -> None: """Add events to template variables.""" for event in self._sbml_events: + trigger = event.getTrigger() + if trigger is None or trigger.getMath() is None: + # An event with no trigger, or a trigger with no MathML, can never fire. + continue self._reject_unsupported_delay(event) label = event.getName().strip() event_type = self._guess_event_type(label) # Compute the trigger formula first: it mutates the trigger AST in place (stripping # units, renaming the avogadro csymbol), and the distance below deep-copies that AST. - trigger_formula = self._formula_to_string(event.getTrigger().getMath()) + trigger_formula = self._formula_to_string(trigger.getMath()) trigger_distance = self._event_trigger_distance(event) assignments = self._event_assignments(event) self._compensate_compartment_resizes(event, assignments) @@ -554,18 +569,31 @@ def _format_events(self) -> None: priority, ) + def _reject_unsupported_packages(self) -> None: + """Raise if the model uses an SBML package the generator cannot translate. + + :raises NotImplementedError: if the model uses the flux-balance-constraints (fbc) package. + """ + if self._sbml_model.getPlugin("fbc") is not None: + raise NotImplementedError("Flux balance constraint models (SBML fbc package) are not supported.") + def _reject_unsupported_delay(self, event) -> None: """Raise if the event has a non-zero delay (delays are not supported).""" - if event.isSetDelay(): - math = self._formula_to_string(event.getDelay().getMath()) - try: - delay = float(math) - except ValueError: - delay = 9999 - - if delay != 0.0: - # Delay of zero is equivalent to no delay - raise NotImplementedError("Events with delays are not supported.") + if not event.isSetDelay(): + return + delay_math = event.getDelay().getMath() + if delay_math is None: + # A delay element with no MathML is treated as no delay. + return + math = self._formula_to_string(delay_math) + try: + delay = float(math) + except ValueError: + delay = 9999 + + if delay != 0.0: + # Delay of zero is equivalent to no delay + raise NotImplementedError("Events with delays are not supported.") @staticmethod def _guess_event_type(label: str) -> EventType: @@ -736,6 +764,9 @@ def _event_priority(self, event) -> Optional[str]: def _format_function_definitions(self) -> None: """Add function definitions to template variables.""" for fd in self._sbml_function_definitions: + if fd.getBody() is None: + # A function definition with no MathML body cannot be emitted. + continue fd_id = fd.getId() label = fd.getName().strip() arg_list = get_function_definition_arguments(fd) @@ -828,8 +859,8 @@ def _format_reactions(self) -> None: # A reaction's ID denotes its rate (flux), an observable that rules, events and # test outputs can read. Expose it as a derived quantity. The flux member is # already declared and computed by the reaction machinery, so it keeps its - # VarType.REACTION type and is not re-declared (is_reaction=True). - self._add_derived_quantity(id_, label, None, NON_DIM_UNITS, is_reaction=True) + # VarType.REACTION type and is not re-declared (kind=REACTION). + self._add_derived_quantity(id_, label, None, NON_DIM_UNITS, kind=DerivedQuantityKind.REACTION) kinetic_law = reaction.getKineticLaw() if kinetic_law is None: @@ -1079,8 +1110,11 @@ def _add_species_dynamics( self._add_equation( var=state_var.derivative_id, math=parseL3Formula(rhs), eq_type=EquationType.DERIVATIVE ) + elif species_id in self._event_assigned_ids: + # Otherwise constant but modified by an event: model as a (variable) parameter + self._add_parameter(species_id, label, initial_value, units) else: - # Truly constant: no reactions, no rules, constant compartment + # Truly constant: no reactions, no rules, no events, constant compartment self._add_derived_quantity(species_id, label, initial_value, units) def _formula_to_string(self, math: "ASTNode", local_parameters: Optional[list["LocalParameter"]] = None) -> str: @@ -1093,7 +1127,7 @@ def _formula_to_string(self, math: "ASTNode", local_parameters: Optional[list["L :param local_parameters: Local parameters in scope for the formula. :return: The equivalent C++ string. """ - return formula_to_string(math, self._variable_types, self._state_variables, local_parameters) + return formula_to_string(math, self._variable_type_map(), self._state_variables, local_parameters) def _get_timescale_multiplier(self) -> float: """Get the timescale multiplier. @@ -1138,17 +1172,38 @@ def _build_variable_index(self) -> dict[str, int]: for record in collection: index_of[record.id] = record.index for derived_quantity in self._derived_quantities: - if not derived_quantity.is_reaction: + if derived_quantity.kind != DerivedQuantityKind.REACTION: index_of[derived_quantity.id] = derived_quantity.index return index_of + def _variable_type_map(self) -> dict[str, VarType]: + """Map each variable id to its VarType, derived from the built record collections. + + The type of a variable is implied by which collection its record belongs to, so no separate + mapping is kept. A reaction flux is both a Reaction (typed VarType.REACTION) and an + REACTION-kind DerivedQuantity; the latter is skipped so the reaction id keeps its type. + """ + types: dict[str, VarType] = {} + for derived_quantity in self._derived_quantities: + if derived_quantity.kind != DerivedQuantityKind.REACTION: + types[derived_quantity.id] = VarType.DERIVED_QUANTITY + for function in self._functions: + types[function.id] = VarType.FUNCTION + for parameter in self._parameters: + types[parameter.id] = VarType.PARAMETER + for reaction in self._reactions: + types[reaction.id] = VarType.REACTION + for state_variable in self._state_variables: + types[state_variable.id] = VarType.STATE_VARIABLE + return types + def _get_variable_type(self, var_id: str) -> VarType: """Get the type of a variable based on its ID. :param var_id: The variable ID. :return: The variable's VarType, or VarType.UNKNOWN if it is not a known variable. """ - return self._variable_types.get(var_id, VarType.UNKNOWN) + return self._variable_type_map().get(var_id, VarType.UNKNOWN) def _order_derived_quantities(self) -> None: """Order derived quantities: normal quantities, then reactions, then conversions. @@ -1161,9 +1216,9 @@ def _order_derived_quantities(self) -> None: """ def group(dq: "DerivedQuantity") -> int: - if dq.is_conversion: + if dq.kind == DerivedQuantityKind.CONVERSION: return 2 - if dq.is_reaction: + if dq.kind == DerivedQuantityKind.REACTION: return 1 return 0 @@ -1173,7 +1228,8 @@ def group(dq: "DerivedQuantity") -> int: def build(self) -> None: """Process the SBML model to set up the formatted variables for templates.""" - self._variable_types = {} + self._reject_unsupported_packages() + self._index_of = None # id -> index cache, built lazily by _get_variable_index self._odes = {} @@ -1191,6 +1247,12 @@ def build(self) -> None: self._events = [] self._functions = [] + # Ids assigned by some event, so species classification can model an otherwise-constant + # event-modified species as a (variable) parameter. + self._event_assigned_ids = { + ea.getVariable() for event in self._sbml_events for ea in event.getListOfEventAssignments() + } + # Names already claimed by real SBML entities and the Chaste base classes. Synthetic # identifiers (derivatives, amount/concentration conversions, initial-assignment # intermediates) are allocated against this so they never collide. @@ -1207,7 +1269,11 @@ def build(self) -> None: self._format_parameters() if not self._state_variables: - raise NotImplementedError("Models without state variables are not supported.") + # The model has no continuous dynamics (no reactions, ODEs or rate rules). Add a + # placeholder state variable with a zero derivative so the generated ODE system still + # has something for the solver to integrate; its outputs (constants, assignment rules + # of time, event-driven changes) are recomputed each step like any other model. + self._add_placeholder_state_variable() self._format_reactions() self._format_function_definitions() @@ -1243,10 +1309,12 @@ def _check_name_conflicts(self) -> None: identifiers.append((var.id, "state variable")) identifiers.append((var.derivative_id, "state-variable derivative")) for dq in self._derived_quantities: - if dq.is_reaction: + if dq.kind == DerivedQuantityKind.REACTION: continue # Declared under reactions; counted there. - kind = "amount/concentration conversion" if dq.is_conversion else "derived quantity" - identifiers.append((dq.id, kind)) + kind_label = ( + "amount/concentration conversion" if dq.kind == DerivedQuantityKind.CONVERSION else "derived quantity" + ) + identifiers.append((dq.id, kind_label)) identifiers += [(s.id, "stoichiometry variable") for s in self._stoichiometry_variables] identifiers += [(r.id, "reaction") for r in self._reactions] identifiers += [(f.id, "function") for f in self._functions] diff --git a/chaste_sbml/_names.py b/chaste_sbml/_names.py index 31635fd8..1c940376 100644 --- a/chaste_sbml/_names.py +++ b/chaste_sbml/_names.py @@ -23,6 +23,8 @@ import re +import libsbml + # C++ keywords (and alternative operator spellings); none may be used as an identifier. CPP_KEYWORDS = frozenset( { @@ -156,8 +158,13 @@ } ) +# C++ macros (from /) that expand when used as a bare identifier, so an SBML id +# spelled like one cannot be emitted verbatim (e.g. a parameter NAN becomes `double NAN;`, which +# the preprocessor rewrites into the macro's expansion). +CPP_MACROS = frozenset({"NAN", "INFINITY", "HUGE_VAL", "EOF", "NULL"}) + # Every name the generator may not use for one of its own identifiers. -RESERVED_NAMES = CPP_KEYWORDS | CHASTE_RESERVED_NAMES +RESERVED_NAMES = CPP_KEYWORDS | CHASTE_RESERVED_NAMES | CPP_MACROS _IDENTIFIER_RE = re.compile(r"^[A-Za-z_][A-Za-z0-9_]*$") @@ -239,7 +246,7 @@ def resolve_cpp_name(base: str, taken) -> str: :return: A C++-safe identifier not present in ``taken``. """ candidate = base - while candidate in CPP_KEYWORDS or candidate in CHASTE_RESERVED_NAMES: + while candidate in CPP_KEYWORDS or candidate in CHASTE_RESERVED_NAMES or candidate in CPP_MACROS: candidate += "_" return unique_name(candidate, taken) @@ -257,6 +264,10 @@ def __init__(self, sbml_model) -> None: """:param sbml_model: The libsbml Model whose ids are being turned into C++.""" self._sbml_model = sbml_model self._taken: set = set() + # C++ identifier (the renamed id) -> original SBML id, for ids renamed to dodge a C++ + # keyword or reserved name. Lets the generated code report a variable under its real SBML + # id even though its C++ identifier had to be escaped. + self._sbml_names: dict[str, str] = {} def resolve_real_id_conflicts(self) -> None: """Rename real SBML ids that are C++ keywords or reserved Chaste names. @@ -289,13 +300,51 @@ def resolve_real_id_conflicts(self) -> None: if not element.isSetId(): continue old_id = element.getId() - if old_id not in CPP_KEYWORDS and old_id not in CHASTE_RESERVED_NAMES: + if old_id not in CPP_KEYWORDS and old_id not in CHASTE_RESERVED_NAMES and old_id not in CPP_MACROS: continue new_id = resolve_cpp_name(old_id, taken) element.setId(new_id) for referrer in model.getListOfAllElements(): referrer.renameSIdRefs(old_id, new_id) + # ``renameSIdRefs`` rewrites variable references but not function-call names (a call + # ``f(x)`` is an ``AST_FUNCTION`` node, not an ``AST_NAME``), so a renamed function + # definition would leave its callers dangling. Rewrite those calls explicitly. + self._rename_function_call_refs(old_id, new_id) taken.add(new_id) + self._sbml_names[new_id] = old_id + + def _rename_function_call_refs(self, old_id: str, new_id: str) -> None: + """Rename every ``old_id(...)`` function call in the model's math to ``new_id(...)``. + + Function calls reference a function definition by name through an ``AST_FUNCTION`` node, + which ``renameSIdRefs`` leaves untouched; without this a function definition renamed to + avoid a keyword/reserved/macro clash (e.g. one called ``this``) keeps being called under + its old, now-invalid name. + """ + + def rename(node: "libsbml.ASTNode") -> None: + if node is None: + return + if node.getType() == libsbml.AST_FUNCTION and node.getName() == old_id: + node.setName(new_id) + for i in range(node.getNumChildren()): + rename(node.getChild(i)) + + for element in self._sbml_model.getListOfAllElements(): + if hasattr(element, "isSetMath") and element.isSetMath(): + rename(element.getMath()) + + def sbml_name(self, cpp_id: str) -> str: + """Return the original SBML id for a (possibly renamed) C++ identifier. + + A variable is reported to Chaste under this name, so it is looked up by its real SBML id + even when the emitted C++ identifier had to be escaped (e.g. a parameter ``time`` emitted + as ``time_``). Ids that were never renamed are returned unchanged. + + :param cpp_id: The C++ identifier as it appears in the generated code. + :return: The original SBML id if ``cpp_id`` was a renamed id, else ``cpp_id``. + """ + return self._sbml_names.get(cpp_id, cpp_id) def reset(self) -> None: """Recompute the taken-name set from the (resolved) model, ready for a build. diff --git a/chaste_sbml/_records.py b/chaste_sbml/_records.py index 1058850b..4abba602 100644 --- a/chaste_sbml/_records.py +++ b/chaste_sbml/_records.py @@ -9,7 +9,7 @@ from dataclasses import dataclass from typing import TYPE_CHECKING, Optional -from ._config import EquationType, EventType, VarType +from ._config import DerivedQuantityKind, EquationType, EventType, VarType if TYPE_CHECKING: from libsbml import ASTNode @@ -45,6 +45,8 @@ class StateVariable: label: str initial_value: Optional[float] units: str + # SBML id reported to Chaste; equals ``id`` unless the id was escaped for C++ (see NameManager). + name: str = "" @dataclass @@ -57,6 +59,8 @@ class Parameter: label: str initial_value: Optional[float] units: str + # SBML id reported to Chaste; equals ``id`` unless the id was escaped for C++ (see NameManager). + name: str = "" @dataclass @@ -67,9 +71,10 @@ class DerivedQuantity: label: str index: int initial_value: Optional[float] - is_conversion: bool - is_reaction: bool units: str + kind: DerivedQuantityKind = DerivedQuantityKind.NORMAL + # SBML id reported to Chaste; equals ``id`` unless the id was escaped for C++ (see NameManager). + name: str = "" @dataclass diff --git a/chaste_sbml/_rendering.py b/chaste_sbml/_rendering.py index fff44b51..29aba91e 100644 --- a/chaste_sbml/_rendering.py +++ b/chaste_sbml/_rendering.py @@ -18,6 +18,7 @@ CONCENTRATION_PREFIX, PREFIX_SEP, ROOT_DIR, + DerivedQuantityKind, EquationType, EventType, VarType, @@ -35,6 +36,7 @@ class CodeRenderer: ) _env.globals["AMOUNT_PREFIX"] = AMOUNT_PREFIX _env.globals["CONCENTRATION_PREFIX"] = CONCENTRATION_PREFIX + _env.globals["DerivedQuantityKind"] = DerivedQuantityKind _env.globals["EquationType"] = EquationType _env.globals["EventType"] = EventType _env.globals["PREFIX_SEP"] = PREFIX_SEP diff --git a/chaste_sbml/templates/cases/semantic.hpp b/chaste_sbml/templates/cases/semantic.hpp index fedd2f0e..36299f24 100644 --- a/chaste_sbml/templates/cases/semantic.hpp +++ b/chaste_sbml/templates/cases/semantic.hpp @@ -160,11 +160,27 @@ class Test{{ model_name }} : public CxxTest::TestSuite TS_ASSERT_EQUALS(values.size(), expected_result_data.size()); for (unsigned i = 0; i < expected_result_data.size(); i++) { - double delta = std::abs(expected_result_data[i][j] - values[i]); - double tol = tol_absolute + tol_relative * std::abs(expected_result_data[i][j]); - std::string msg(sth::ToString(values[i]) + " vs " + sth::ToString(expected_result_data[i][j]) + double expected = expected_result_data[i][j]; + double actual = values[i]; + std::string msg(sth::ToString(actual) + " vs " + sth::ToString(expected) + " at " + sth::ToString(ode_solution.rGetTimes()[i], 3) + " for " + var_name); - TSM_ASSERT_LESS_THAN_EQUALS(msg.c_str(), delta, tol); + // Non-finite expected values (e.g. a division by zero gives +/-INF, 0/0 gives + // NaN) cannot be checked with a tolerance: |INF - INF| is NaN, and any + // comparison against NaN is false. Match them by kind (and sign) instead. + if (std::isnan(expected)) + { + TSM_ASSERT(msg.c_str(), std::isnan(actual)); + } + else if (std::isinf(expected)) + { + TSM_ASSERT(msg.c_str(), std::isinf(actual) && std::signbit(actual) == std::signbit(expected)); + } + else + { + double delta = std::abs(expected - actual); + double tol = tol_absolute + tol_relative * std::abs(expected); + TSM_ASSERT_LESS_THAN_EQUALS(msg.c_str(), delta, tol); + } } } diff --git a/chaste_sbml/templates/ode/ode.cpp b/chaste_sbml/templates/ode/ode.cpp index e22ee65d..93aadb4d 100644 --- a/chaste_sbml/templates/ode/ode.cpp +++ b/chaste_sbml/templates/ode/ode.cpp @@ -327,7 +327,7 @@ void CellwiseOdeSystemInformation<{{ ode_class_name }}>::Initialise() { // STATE VARIABLES {% for var in state_variables %} - this->mVariableNames.push_back("{{ var['id'] }}"); + this->mVariableNames.push_back("{{ var['name'] }}"); this->mVariableUnits.push_back("{{ var['units'] }}"); {% if var["initial_value"] is not none %} this->mInitialConditions.push_back({{ var['initial_value'] }}); @@ -339,14 +339,14 @@ void CellwiseOdeSystemInformation<{{ ode_class_name }}>::Initialise() // DERIVED QUANTITIES {% for dq in derived_quantities %} - this->mDerivedQuantityNames.push_back("{{ dq['id'] }}"); + this->mDerivedQuantityNames.push_back("{{ dq['name'] }}"); this->mDerivedQuantityUnits.push_back("{{ dq['units'] }}"); {% endfor %} // PARAMETERS {% for var in parameters %} - this->mParameterNames.push_back("{{ var['id'] }}"); + this->mParameterNames.push_back("{{ var['name'] }}"); this->mParameterUnits.push_back("{{ var['units'] }}"); {% endfor %} diff --git a/chaste_sbml/templates/ode/ode.hpp b/chaste_sbml/templates/ode/ode.hpp index 76d03632..a84752ae 100644 --- a/chaste_sbml/templates/ode/ode.hpp +++ b/chaste_sbml/templates/ode/ode.hpp @@ -42,7 +42,7 @@ class {{ ode_class_name }} : public AbstractSbmlOdeSystem // DERIVED QUANTITIES {% for dq in derived_quantities %} -{% if dq["is_conversion"] is false() and dq["is_reaction"] is false() %} +{% if dq["kind"] == DerivedQuantityKind.NORMAL %} double {{ dq["id"] }}; // {{ dq["label"] }} {% endif %} {% endfor %} diff --git a/chaste_sbml/tests/test_expressions.py b/chaste_sbml/tests/test_expressions.py index 24ad0fbd..70d57838 100644 --- a/chaste_sbml/tests/test_expressions.py +++ b/chaste_sbml/tests/test_expressions.py @@ -1,5 +1,6 @@ """Tests for SBML expression/AST -> C++ formula translation.""" +import libsbml import pytest from libsbml import ( AST_FUNCTION_COS, @@ -15,10 +16,12 @@ from chaste_sbml._config import VarType from chaste_sbml._expressions import ( + AVOGADRO_PLACEHOLDER, collect_ast_names, - convert_infix_operator_to_function_syntax, formula_to_string, - replace_avogadro_csymbol, + replace_constant_symbols, + rewrite_nary_relational, + rewrite_power, search_ast_type, strip_ast_units, substitute_ast_names, @@ -26,44 +29,43 @@ from chaste_sbml._records import StateVariable -def test_convert_infix_operator_to_function_syntax_power(): - """Converts infix power expressions to pow calls.""" - converted = convert_infix_operator_to_function_syntax( - formula="(a + b) ^ (1.0 / 2.0)", - operator="^", - function_name="pow", - ) - - assert converted == "pow((a + b), (1.0 / 2.0))" - - -def test_convert_infix_operator_to_function_syntax_custom_pair(): - """Supports arbitrary operator/function conversion pairs.""" - converted = convert_infix_operator_to_function_syntax( - formula="alpha @@ beta", operator="@@", function_name="combine" - ) +@pytest.mark.parametrize( + ("formula", "expected"), + [ + ("a ^ b", "std::pow(a, b)"), + ("2 ^ 3", "std::pow(2.0, 3.0)"), + ("(a + b) ^ 2", "std::pow(a + b, 2.0)"), + # Nested power fully converts (the old string rewrite left an inner ^ behind). + ("a ^ b ^ c", "std::pow(a, std::pow(b, c))"), + ], +) +def test_formula_to_string_power_uses_std_pow(formula, expected): + """Power operators (including nested) become std::pow calls.""" + variable_types = {name: VarType.PARAMETER for name in ("a", "b", "c")} + assert formula_to_string(parseL3Formula(formula), variable_types, state_variables=[]) == expected - assert converted == "combine(alpha, beta)" +def test_rewrite_power_converts_only_power_nodes(): + """rewrite_power converts AST_POWER nodes to pow() and leaves other nodes alone.""" + node = parseL3Formula("a ^ b ^ c") + rewrite_power(node) + assert formulaToL3String(node) == "pow(a, pow(b, c))" -def test_convert_infix_operator_to_function_syntax_nested_and_repeated(): - """Converts nested and repeated operator expressions in one formula.""" - converted = convert_infix_operator_to_function_syntax(formula="(x ^ y) ^ z", operator="^", function_name="pow") + unaffected = parseL3Formula("a * b + c") + rewrite_power(unaffected) + assert formulaToL3String(unaffected) == "a * b + c" - assert converted == "pow((pow(x, y)), z)" +def test_rewrite_power_converts_mathml_function_power(): + """A MathML (AST_FUNCTION_POWER, distinct from L3 AST_POWER) is also converted. -@pytest.mark.parametrize( - ("operator", "function_name", "error_message"), - [ - ("", "pow", "operator must be a non-empty string"), - ("^", "", "function_name must be a non-empty string"), - ], -) -def test_convert_infix_operator_to_function_syntax_invalid_args(operator: str, function_name: str, error_message: str): - """Validates required arguments for operator/function conversion.""" - with pytest.raises(ValueError, match=error_message): - convert_infix_operator_to_function_syntax(formula="a ^ b", operator=operator, function_name=function_name) + Guards against only handling the L3-parsed power node type: SBML MathML powers parse to a + different AST type that ``formulaToL3String`` still renders as ``^``. + """ + mathml = '' "a2" + node = libsbml.readMathMLFromString(mathml) + rewrite_power(node) + assert formulaToL3String(node) == "pow(a, 2)" def test_formula_to_string_functions_constants_and_power(): @@ -75,6 +77,42 @@ def test_formula_to_string_functions_constants_and_power(): assert out == "k * std::pow(S, 2.0) + M_PI" +def test_formula_to_string_implies_maps_to_sbmlmath(): + """The MathML 'implies' logical operator maps to the sm::implies helper.""" + assert formula_to_string(parseL3Formula("implies(a, b)"), {}, state_variables=[]) == "sm::implies(a, b)" + + +@pytest.mark.parametrize( + ("formula", "expected"), + [ + # n-ary relationals (>2 operands) become sm:: calls with the correct semantics... + ("lt(1, 2, 3)", "sm::lt(1.0, 2.0, 3.0)"), + ("gt(2, 1, 2)", "sm::gt(2.0, 1.0, 2.0)"), + ("leq(1, 2, 3, 4)", "sm::leq(1.0, 2.0, 3.0, 4.0)"), + ("geq(2, 1, 2)", "sm::geq(2.0, 1.0, 2.0)"), + ("eq(1, 1, 2)", "sm::eq(1.0, 1.0, 2.0)"), + # ...while binary relationals stay as infix operators (unchanged output). + ("a < b", "a < b"), + ("a >= b", "a >= b"), + ], +) +def test_formula_to_string_nary_relational_uses_sbmlmath(formula, expected): + """n-ary relationals route to the variadic sm:: helpers; binary ones stay infix.""" + variable_types = {"a": VarType.PARAMETER, "b": VarType.PARAMETER} + assert formula_to_string(parseL3Formula(formula), variable_types, state_variables=[]) == expected + + +def test_rewrite_nary_relational_leaves_binary_untouched(): + """rewrite_nary_relational converts only >2-operand relationals, leaving binary ones as-is.""" + binary = parseL3Formula("a < b") + rewrite_nary_relational(binary) + assert formulaToL3String(binary) == "a < b" + + nary = parseL3Formula("lt(a, b, c)") + rewrite_nary_relational(nary) + assert formulaToL3String(nary) == "lt(a, b, c)" + + def test_formula_to_string_rateof_state_variable_uses_derivative(): """rateOf(state var) resolves to that variable's derivative id.""" math = parseL3Formula("rateOf(X)") @@ -121,8 +159,8 @@ def test_substitute_ast_names_leaves_original_unchanged(): assert formulaToL3String(original) == "a + b" -def test_replace_avogadro_csymbol_renames_nested_node(): - """An avogadro csymbol anywhere in the tree becomes a plain name with the placeholder id.""" +def test_replace_constant_symbols_renames_nested_node(): + """An avogadro csymbol anywhere in the tree becomes a plain name with its placeholder id.""" one = ASTNode(AST_REAL) one.setValue(1.0) avogadro = ASTNode(AST_NAME_AVOGADRO) @@ -130,21 +168,83 @@ def test_replace_avogadro_csymbol_renames_nested_node(): expr.addChild(one) expr.addChild(avogadro) - replace_avogadro_csymbol(expr, "AVO") + replace_constant_symbols(expr) renamed = expr.getChild(1) assert renamed.getType() == AST_NAME - assert renamed.getName() == "AVO" + assert renamed.getName() == AVOGADRO_PLACEHOLDER -def test_replace_avogadro_csymbol_leaves_other_names(): - """A non-avogadro name node is left untouched.""" +def test_replace_constant_symbols_leaves_other_names(): + """A non-symbol name node is left untouched.""" name = ASTNode(AST_NAME) name.setName("x") - replace_avogadro_csymbol(name, "AVO") + replace_constant_symbols(name) assert name.getName() == "x" +@pytest.mark.parametrize( + ("mathml_symbol", "expected"), + [ + ("", "M_PI"), + ("", "M_E"), + ("", "std::numeric_limits::infinity()"), + ("", "std::numeric_limits::quiet_NaN()"), + ], +) +def test_formula_to_string_math_constants(mathml_symbol, expected): + """SBML math constants map to their C++ form.""" + math = libsbml.readMathMLFromString(f'{mathml_symbol}') + assert formula_to_string(math, {}, state_variables=[]) == expected + + +@pytest.mark.parametrize( + ("mathml_symbol", "param", "expected"), + [ + ("", "INF", "std::numeric_limits::infinity()"), + ("", "NaN", "std::numeric_limits::quiet_NaN()"), + ], +) +def test_formula_to_string_infinity_nan_distinct_from_parameter(mathml_symbol, param, expected): + """The / constants stay distinct from a parameter spelled INF/NaN. + + libsbml stores these constants as real nodes that render as ``INF``/``NaN`` -- the same text as + a reference to a parameter so named -- so the constant must be resolved from the AST, not its + rendered spelling (test-suite cases 1811/1813). + """ + math = libsbml.readMathMLFromString(f'{mathml_symbol}') + assert formula_to_string(math, {param: VarType.PARAMETER}, state_variables=[]) == expected + + +def test_formula_to_string_pi_constant_and_parameter_are_distinct(): + """The constant maps to M_PI while a same-named parameter keeps its own name.""" + math = libsbml.readMathMLFromString( + 'pi' + ) + assert formula_to_string(math, {"pi": VarType.PARAMETER}, state_variables=[]) == "pi + M_PI" + + +def test_formula_to_string_time_symbol_and_named_variable_are_distinct(): + """The time csymbol maps to `time`; a plain variable named `t` is left as-is (not time).""" + assert formula_to_string(parseL3Formula("time"), {}, state_variables=[]) == "time" + assert formula_to_string(parseL3Formula("t"), {"t": VarType.PARAMETER}, state_variables=[]) == "t" + + +def test_formula_to_string_rateof_local_parameter_is_zero(): + """rateOf of a kinetic-law local parameter (constant) is zero, not a global's derivative.""" + from chaste_sbml._records import LocalParameter + + out = formula_to_string( + parseL3Formula("rateOf(k)"), + variable_types={"k": VarType.STATE_VARIABLE}, + state_variables=[ + StateVariable(index=0, id="k", derivative_id="d_k_dt", label="", initial_value=None, units="") + ], + local_parameters=[LocalParameter(id="k", label="", value="1")], + ) + assert out == "0.0" + + def test_strip_ast_units_removes_units_recursively(): """Units annotations are stripped from a node and its children.""" node = parseL3Formula("2 * 3") diff --git a/chaste_sbml/tests/test_model_builder.py b/chaste_sbml/tests/test_model_builder.py index 8822ef87..ea76aff2 100644 --- a/chaste_sbml/tests/test_model_builder.py +++ b/chaste_sbml/tests/test_model_builder.py @@ -4,7 +4,15 @@ import pytest from libsbml import parseL3Formula -from chaste_sbml._config import ROOT_DIR, EquationType, VarType +from chaste_sbml._config import ( + CHASTE_PREFIX, + PLACEHOLDER_STATE_ID, + PREFIX_SEP, + ROOT_DIR, + DerivedQuantityKind, + EquationType, + VarType, +) from chaste_sbml._model_builder import ModelBuilder from chaste_sbml._names import NameConflictError, NameManager from chaste_sbml._records import ( @@ -94,10 +102,8 @@ def _make_state_variable(*, id, derivative_id): return StateVariable(index=0, id=id, derivative_id=derivative_id, label="", initial_value=None, units="") -def _make_derived_quantity(*, id, is_conversion=False, is_reaction=False): - return DerivedQuantity( - id=id, label="", index=0, initial_value=None, is_conversion=is_conversion, is_reaction=is_reaction, units="" - ) +def _make_derived_quantity(*, id, kind=DerivedQuantityKind.NORMAL): + return DerivedQuantity(id=id, label="", index=0, initial_value=None, units="", kind=kind) def test_check_name_conflicts_passes_for_clean_model(): @@ -110,7 +116,7 @@ def test_check_name_conflicts_flags_synthetic_vs_real_collision(): """A real id equal to a synthesised amount-conversion name is caught.""" builder = _builder_with_names( state_variables=["X"], - derived_quantities=[{"id": "amt__X", "is_conversion": True, "is_reaction": False}], + derived_quantities=[{"id": "amt__X", "kind": DerivedQuantityKind.CONVERSION}], parameters=["amt__X"], ) with pytest.raises(NameConflictError, match="amt__X"): @@ -136,7 +142,7 @@ def test_check_name_conflicts_ignores_reaction_flux_duplicate(): builder = _builder_with_names( state_variables=["C"], reactions=["reaction1"], - derived_quantities=[{"id": "reaction1", "is_conversion": False, "is_reaction": True}], + derived_quantities=[{"id": "reaction1", "kind": DerivedQuantityKind.REACTION}], ) builder._check_name_conflicts() # should not raise @@ -171,6 +177,118 @@ def test_check_name_conflicts_flags_local_parameter_shadowing_time(): builder._check_name_conflicts() +PLACEHOLDER_ID = PREFIX_SEP.join([CHASTE_PREFIX, PLACEHOLDER_STATE_ID]) + + +def _no_ode_model(tmp_path): + """Write a no-ODE SBML model (constant species + time-dependent assignment rule) to disk. + + The model has no reactions, ODEs or rate rules, so it has no state variables of its own. + """ + doc = libsbml.SBMLDocument(3, 2) + model = doc.createModel() + model.setId("NoOde") + compartment = model.createCompartment() + compartment.setId("cell") + compartment.setConstant(True) + compartment.setSize(1.0) + compartment.setSpatialDimensions(3) + # A parameter that is an explicit function of time via an assignment rule. + parameter = model.createParameter() + parameter.setId("P") + parameter.setConstant(False) + parameter.setValue(0.0) + rule = model.createAssignmentRule() + rule.setVariable("P") + rule.setMath(parseL3Formula("2 * time")) + # A constant species. + species = model.createSpecies() + species.setId("S") + species.setCompartment("cell") + species.setConstant(True) + species.setBoundaryCondition(True) + species.setHasOnlySubstanceUnits(False) + species.setInitialConcentration(3.0) + + path = tmp_path / "NoOde.xml" + libsbml.writeSBMLToFile(doc, str(path)) + return path + + +def test_build_synthesises_placeholder_state_variable_for_no_ode_model(tmp_path): + """A model with no ODEs gets a single zero-derivative placeholder state variable.""" + builder = _build(_no_ode_model(tmp_path)) + + assert [s.id for s in builder._state_variables] == [PLACEHOLDER_ID] + placeholder = builder._state_variables[0] + assert placeholder.initial_value == 0.0 + + deriv_equations = [e for e in builder._equations if e.type == EquationType.DERIVATIVE] + assert len(deriv_equations) == 1 + assert deriv_equations[0].var == placeholder.derivative_id + assert deriv_equations[0].rhs == "0.0" + + # The placeholder's member must be initialised to zero: Initialise() reads it to set the + # solver's initial condition, so without this equation the state would be uninitialised memory. + init_equations = [e for e in builder._equations if e.type == EquationType.INITIAL_VALUE] + assert placeholder.id in [e.var for e in init_equations] + assert next(e for e in init_equations if e.var == placeholder.id).rhs == "0.0" + + +def test_build_no_ode_model_preserves_real_outputs(tmp_path): + """The placeholder does not displace the model's real outputs (assignment rule, species).""" + builder = _build(_no_ode_model(tmp_path)) + + # The time-dependent assignment-rule parameter is still emitted as a derived quantity, and + # the constant boundary species is still emitted (as a parameter). Neither is the placeholder. + assert "P" in {d.id for d in builder._derived_quantities} + assert "S" in {p.id for p in builder._parameters} + + +def test_build_event_modified_constant_species_is_a_parameter(tmp_path): + """An otherwise-constant species modified by an event is a parameter, not a derived quantity. + + A derived quantity is recomputed from the current member each step, so an event's change would + be reported at every earlier time point; a parameter is time-resolved in the recorded solution. + """ + doc = libsbml.SBMLDocument(3, 2) + model = doc.createModel() + compartment = model.createCompartment() + compartment.setId("C") + compartment.setConstant(True) + compartment.setSize(1.0) + compartment.setSpatialDimensions(3) + species = model.createSpecies() + species.setId("S") + species.setCompartment("C") + species.setConstant(False) + species.setBoundaryCondition(False) + species.setHasOnlySubstanceUnits(False) + species.setInitialAmount(0.0) + event = model.createEvent() + event.setUseValuesFromTriggerTime(True) + trigger = event.createTrigger() + trigger.setMath(parseL3Formula("time >= 1")) + trigger.setInitialValue(True) + trigger.setPersistent(True) + assignment = event.createEventAssignment() + assignment.setVariable("S") + assignment.setMath(parseL3Formula("5")) + path = tmp_path / "event_species.xml" + libsbml.writeSBMLToFile(doc, str(path)) + + builder = _build(path) + + assert "S" in {p.id for p in builder._parameters} + assert "S" not in {d.id for d in builder._derived_quantities} + + +def test_build_does_not_add_placeholder_when_model_has_odes(): + """A model with genuine ODEs is left untouched: no placeholder state variable is added.""" + builder = _build(REFERENCE / "Goldbeter1991" / "Goldbeter1991.xml") + assert PLACEHOLDER_ID not in {s.id for s in builder._state_variables} + + def test_build_extracts_state_variables_and_derivatives(): """Each SBML species with an ODE becomes a state variable with a d__dt derivative.""" builder = _build(REFERENCE / "Goldbeter1991" / "Goldbeter1991.xml") @@ -188,7 +306,7 @@ def test_build_extracts_parameters_and_reactions(): def test_build_adds_amount_conversions_for_amount_species(): """Amount species get an amt__ conversion derived quantity.""" builder = _build(REFERENCE / "Goldbeter1991" / "Goldbeter1991.xml") - conversions = {d.id for d in builder._derived_quantities if d.is_conversion} + conversions = {d.id for d in builder._derived_quantities if d.kind == DerivedQuantityKind.CONVERSION} assert conversions == {"amt__C", "amt__M", "amt__X"} @@ -208,14 +326,96 @@ def test_build_extracts_cell_division_event(): assert event.assignments == [EventAssignment(index=3, lhs="m", rhs="m / 2.0", type=VarType.STATE_VARIABLE)] +def _event_builder(): + """A bare builder with an empty model, wired just enough for _format_events.""" + builder = _builder_without_init() + doc = libsbml.SBMLDocument(3, 2) + builder._sbml_model = doc.createModel() + builder._keep_doc = doc # keep the owning document (and thus the model) alive + builder._sbml_events = builder._sbml_model.getListOfEvents() + builder._sbml_function_definitions = builder._sbml_model.getListOfFunctionDefinitions() + builder._state_variables = [] + builder._parameters = [] + builder._derived_quantities = [] + builder._reactions = [] + builder._events = [] + builder._functions = [] + return builder + + +def test_format_events_skips_event_with_missing_trigger_math(): + """An event whose trigger has no MathML can never fire, so no event is emitted.""" + builder = _event_builder() + event = builder._sbml_model.createEvent() + event.createTrigger() # trigger element present but with no math + event.createEventAssignment().setVariable("p") + + builder._format_events() + + assert builder._events == [] + + +def test_format_events_skips_event_with_no_trigger(): + """An event with no trigger element at all can never fire, so no event is emitted.""" + builder = _event_builder() + builder._sbml_model.createEvent() # no trigger + + builder._format_events() + + assert builder._events == [] + + +def test_reject_unsupported_delay_allows_missing_delay_math(): + """A delay element with no MathML is treated as no delay, so it is not rejected.""" + builder = _event_builder() + event = builder._sbml_model.createEvent() + event.createDelay() # delay element present but with no math + + builder._reject_unsupported_delay(event) # should not raise + + +def test_reject_unsupported_packages_rejects_fbc(): + """A model using the flux-balance-constraints (fbc) package is rejected.""" + doc = libsbml.SBMLDocument(libsbml.SBMLNamespaces(3, 1, "fbc", 2)) + builder = _builder_without_init() + builder._sbml_model = doc.createModel() + builder._keep_doc = doc + + with pytest.raises(NotImplementedError, match="fbc"): + builder._reject_unsupported_packages() + + +def test_reject_unsupported_packages_allows_plain_model(): + """A model without an unsupported package is not rejected.""" + doc = libsbml.SBMLDocument(3, 2) + builder = _builder_without_init() + builder._sbml_model = doc.createModel() + builder._keep_doc = doc + + builder._reject_unsupported_packages() # should not raise + + +def test_format_function_definitions_skips_missing_body(): + """A function definition with no MathML body is skipped rather than crashing.""" + builder = _event_builder() + builder._sbml_model.createFunctionDefinition().setId("foo") # no body + + builder._format_function_definitions() + + assert builder._functions == [] + + def _distance_builder() -> ModelBuilder: """A bare builder wired just enough for _event_trigger_distance and trigger formatting.""" builder = _builder_without_init() doc = libsbml.SBMLDocument(3, 2) builder._sbml_model = doc.createModel() builder._keep_doc = doc # keep the owning document (and thus the model) alive - builder._variable_types = {} builder._state_variables = [] + builder._parameters = [] + builder._derived_quantities = [] + builder._reactions = [] + builder._functions = [] return builder @@ -266,8 +466,11 @@ def _resize_builder(*, assignment_rules=()) -> ModelBuilder: builder._sbml_compartments = model.getListOfCompartments() builder._sbml_species = model.getListOfSpecies() builder._assignment_rules = list(assignment_rules) - builder._variable_types = {"S": VarType.STATE_VARIABLE} - builder._state_variables = [] + builder._state_variables = [_make_state_variable(id="S", derivative_id="d_S_dt")] + builder._parameters = [] + builder._derived_quantities = [] + builder._reactions = [] + builder._functions = [] builder._index_of = {"S": 0} return builder diff --git a/chaste_sbml/tests/test_names.py b/chaste_sbml/tests/test_names.py index a188e504..ca4bbf28 100644 --- a/chaste_sbml/tests/test_names.py +++ b/chaste_sbml/tests/test_names.py @@ -166,6 +166,72 @@ def test_name_manager_resolve_leaves_safe_ids_untouched(): assert sbml_model.getElementBySId("cell") is not None +def test_name_manager_sbml_name_recovers_original_id(): + """After a keyword id is renamed, sbml_name maps the C++ id back to the original SBML id. + + A variable must be reported to Chaste under its real SBML id even though the emitted C++ + identifier had to be escaped, so it can be looked up by that id. + """ + doc = libsbml.SBMLDocument(3, 2) + sbml_model = doc.createModel() + parameter = sbml_model.createParameter() + parameter.setId("true") # 'true' is a C++ keyword, so it is renamed to 'true_' + parameter.setConstant(True) + parameter.setValue(1.0) + + manager = NameManager(sbml_model) + manager.resolve_real_id_conflicts() + + assert sbml_model.getParameter("true_") is not None + assert manager.sbml_name("true_") == "true" # renamed id maps back to the real SBML id + assert manager.sbml_name("k") == "k" # an un-renamed id is returned unchanged + + +def test_name_manager_resolve_renames_cpp_macro(): + """An id spelling a C++ macro (e.g. NAN) is renamed but still reported under its original id. + + ``double NAN;`` would be rewritten by the preprocessor into the macro's expansion, so the id + must be escaped for the emitted code while keeping ``NAN`` as its reported name. + """ + doc = libsbml.SBMLDocument(3, 2) + sbml_model = doc.createModel() + parameter = sbml_model.createParameter() + parameter.setId("NAN") # the NAN macro from + parameter.setConstant(True) + parameter.setValue(0.1) + + manager = NameManager(sbml_model) + manager.resolve_real_id_conflicts() + + assert sbml_model.getParameter("NAN") is None + assert sbml_model.getParameter("NAN_") is not None + assert manager.sbml_name("NAN_") == "NAN" + + +def test_name_manager_resolve_renames_function_calls(): + """A function definition renamed to dodge a keyword also has its callers rewritten. + + A call ``this(x)`` is an ``AST_FUNCTION`` node, which ``renameSIdRefs`` leaves untouched, so + renaming the ``this`` function definition must explicitly rewrite the call inside ``caller`` -- + otherwise the emitted C++ would call a function that no longer exists (and ``this`` is a keyword). + """ + doc = libsbml.SBMLDocument(3, 2) + sbml_model = doc.createModel() + callee = sbml_model.createFunctionDefinition() + callee.setId("this") # 'this' is a C++ keyword + callee.setMath(libsbml.parseL3Formula("lambda(x, x)")) + caller = sbml_model.createFunctionDefinition() + caller.setId("caller") + caller.setMath(libsbml.parseL3Formula("lambda(x, this(x))")) + + NameManager(sbml_model).resolve_real_id_conflicts() + + assert sbml_model.getFunctionDefinition("this") is None + assert sbml_model.getFunctionDefinition("this_") is not None + # The body of 'caller' now calls the renamed function. + assert "this_(x)" in libsbml.formulaToL3String(sbml_model.getFunctionDefinition("caller").getMath()) + + def test_generate_header_guard(): """Test header guard generation.""" test_cases = [