diff --git a/Project.toml b/Project.toml index 500a41047..bea7f6502 100644 --- a/Project.toml +++ b/Project.toml @@ -43,7 +43,7 @@ Logging = "1.10, 1.11" MacroTools = "0.5" ModelingToolkitBase = "1.12" Nemo = "0.46, 0.47, 0.48, 0.49, 0.50, 0.51, 0.52, 0.53, 0.54, 0.55, 0.56" -ParamPunPam = "0.5.5" +ParamPunPam = "0.5" PrecompileTools = "1.2" Primes = "0.5" Random = "1.10, 1.11" diff --git a/benchmarking/IdentifiableFunctions/.gitignore b/benchmarking/IdentifiableFunctions/.gitignore index fbca22537..cd07ff2f3 100644 --- a/benchmarking/IdentifiableFunctions/.gitignore +++ b/benchmarking/IdentifiableFunctions/.gitignore @@ -1 +1,2 @@ results/ +logs/ diff --git a/benchmarking/IdentifiableFunctions/benchmark_result.md b/benchmarking/IdentifiableFunctions/benchmark_result.md index 9e46740c0..9a8cdbc9f 100644 --- a/benchmarking/IdentifiableFunctions/benchmark_result.md +++ b/benchmarking/IdentifiableFunctions/benchmark_result.md @@ -1,84 +1,92 @@ ## Benchmark results -2023-09-15T02:10:37.348 +2026-07-28T01:16:05.923 - - Benchmarked function: `find_identifiable_functions` - - Workers: 7 - - Timeout: 600 s +- Benchmarked function: `find_identifiable_functions` +- Workers: 30 +- Timeout: 7200 s +- Memory limit: 28192 MiB -**All timings in seconds.** - -| Model | (:normalforms, 2)_with_states / Runtime | (:normalforms, 2)_with_states / beautifulization | -|:----------------------------- |:--------------------------------------- |:------------------------------------------------ | -| Akt pathway | 2.44 | 0.28 | -| Bilirubin2_io | 4.08 | 0.76 | -| Biohydrogenation_io | 0.44 | 0.11 | -| Bruno2016 | 0.19 | 0.07 | -| CD8 T cell differentiation | 0.72 | 0.27 | -| CGV1990 | 2.25 | 0.31 | -| Chemical reaction network | 0.36 | 0.13 | -| Crauste_SI | 0.47 | 0.18 | -| Fujita | 2.50 | 0.30 | -| Goodwin oscillator | 0.29 | 0.09 | -| HIV | 0.78 | 0.19 | -| HIV2_io | 2.00 | 0.27 | -| HighDimNonLin | 20.55 | 0.68 | -| KD1999 | 0.97 | 0.14 | -| LLW1987_io | 0.25 | 0.08 | -| MAPK model (5 outputs) | 48.45 | 0.37 | -| MAPK model (6 outputs) | 7.99 | 0.42 | -| Modified LV for testing | 0.07 | 0.01 | -| PK1 | 0.74 | 0.18 | -| PK2 | 130.96 | 0.07 | -| Pharm | 152.26 | 0.09 | -| Phosphorylation | 0.45 | 0.11 | -| Pivastatin | 8.44 | 0.04 | -| QY | 49.21 | 0.23 | -| Ruminal lipolysis | 0.21 | 0.07 | -| SEAIJRC Covid model | 85.13 | 0.11 | -| SEIR 34 | 0.39 | 0.09 | -| SEIR 36 ref | 0.78 | 0.14 | -| SEIR2T | 0.23 | 0.09 | -| SEIRT | 0.36 | 0.08 | -| SEIR_1_io | 0.56 | 0.10 | -| SEUIR | - | - | -| SIR 19 | - | - | -| SIR 21 | 0.30 | 0.19 | -| SIR 24 | 0.36 | 0.07 | -| SIR 6 | 0.20 | 0.08 | -| SIRS forced | 10.62 | 0.40 | -| SIWR original | 13.72 | 0.07 | -| SIWR with extra output | 0.67 | 0.10 | -| SLIQR | 1.42 | 0.15 | -| St | 37.54 | 1.38 | -| Transfection_4State | 0.23 | 0.07 | -| Treatment_io | 0.44 | 0.14 | -| cLV1 (2o) | 0.80 | 0.20 | -| cholera | 0.66 | 0.09 | -| generalizedLoktaVolterra (1o) | 0.18 | 0.05 | -| generalizedLoktaVolterra (2o) | 0.09 | 0.07 | -| p53 | 1.21 | 0.34 | +|Model|funcs (cmp_default)|funcs (cmp_lie)| +|---|---|---| +|Akt pathway|[reaction_8_k1, reaction_7_k1, reaction_6_k1, reaction_5_k2, reaction_4_k1, reaction_3_k1, reaction_2_k2, reaction_1_k1 - reaction_1_k2 - reaction_9_k1, reaction_2_k1//reaction_5_k1, a3//reaction_5_k1, a2//reaction_5_k1, a1//reaction_5_k1, pS6(t)*reaction_5_k1, pAkt_S6(t)*reaction_5_k1, S6(t)*reaction_5_k1, pAkt(t)*reaction_5_k1, Akt(t)*reaction_5_k1, pEGFR_Akt(t)*reaction_5_k1, pEGFR(t)*reaction_5_k1, EGF_EGFR(t)*reaction_5_k1*reaction_9_k1]|[reaction_8_k1, reaction_7_k1, reaction_6_k1, reaction_5_k2, reaction_4_k1, reaction_3_k1, reaction_2_k2, reaction_1_k1 - reaction_1_k2 - reaction_9_k1, reaction_2_k1//reaction_5_k1, a3//reaction_5_k1, a2//reaction_5_k1, a1//reaction_5_k1, pS6(t)*reaction_5_k1, pAkt(t)*a2 + pAkt_S6(t)*a2, pEGFR(t)*a1 + pEGFR_Akt(t)*a1, EGF_EGFR(t)*reaction_5_k1*reaction_9_k1, S6(t)*reaction_5_k1, pAkt_S6(t)*reaction_5_k1, Akt(t)*reaction_5_k1, pEGFR_Akt(t)*reaction_5_k1]| +|Bilirubin2_io|[k01, k21 + k31 + k41, k12 + k13 + k14, k21*k31 + k21*k41 + k31*k41, k12*k13 + k12*k14 + k13*k14, k12*k31 + k12*k41 + k13*k21 + k13*k41 + k14*k21 + k14*k31, k21*k31*k41, k12*k13*k14, x1(t), x2(t) + x3(t) + x4(t), x2(t)*x3(t) + x2(t)*x4(t) + x3(t)*x4(t), x2(t)*k31 + x2(t)*k41 + x3(t)*k21 + x3(t)*k41 + x4(t)*k21 + x4(t)*k31, x2(t)*k13 + x2(t)*k14 + x3(t)*k12 + x3(t)*k14 + x4(t)*k12 + x4(t)*k13]|[k01, k21 + k31 + k41, k12 + k13 + k14, k21*k31 + k21*k41 + k31*k41, k12*k13 + k12*k14 + k13*k14, k12*k31 + k12*k41 + k13*k21 + k13*k41 + k14*k21 + k14*k31, k21*k31*k41, k12*k13*k14, x2(t) + x3(t) + x4(t), x1(t), x2(t)*k31 + x2(t)*k41 + x3(t)*k21 + x3(t)*k41 + x4(t)*k21 + x4(t)*k31, x2(t)*k13 + x2(t)*k14 + x3(t)*k12 + x3(t)*k14 + x4(t)*k12 + x4(t)*k13]| +|Biohydrogenation_io|[k7, k6, k5, k10 + 2*k8, k9^2, k10*k9, x5(t), x4(t), x6(t) + k8]|[k7, k6, k5, k10 + 2*k8, k9^2, k10*k9, x4(t), x5(t), x6(t) - 1//2*k10]| +|Bruno2016|[kbeta10, kbeta, kcryOH + kcrybeta, beta10(t), beta(t), cry(t)*kcryOH]|[kbeta10, kbeta, kcryOH + kcrybeta, beta(t), beta10(t), cry(t)*kcryOH]| +|CD8 T cell differentiation|[mu_PL, mu_PE, mu_P, mu_N, mu_M, mu_LL, mu_LE, mu_EE, delta_LM, delta_EL, rho_E//rho_P, delta_NE//rho_P, M(t), S(t), E(t), N(t), P(t)*rho_P]|[mu_PL, mu_PE, mu_P, mu_N, mu_M, mu_LL, mu_LE, mu_EE, delta_LM, delta_EL, rho_E//rho_P, delta_NE//rho_P, M(t), N(t), S(t), E(t), P(t)*rho_P]| +|CGV1990|[k7, k6, k4, k3, R*V36 + 1//25*S*V3, R*V3 + S*V36, R*V36*k5 + 1//5*S*V36*k5, (V3*k5 + 5*V36*k5)//V3, q7(t), q3(t), q1(t), q35(t) + q36(t), q35(t)*q36(t), q35(t)*S*V36 + q36(t)*R*V3]|[k7, k6, k4, k3, R*V36 + 1//25*S*V3, R*V3 + S*V36, R*V36*k5 + 1//5*S*V36*k5, (V3*k5 + 5*V36*k5)//V3, q7(t), q1(t), q35(t)*S*V36 + q36(t)*R*V3, q35(t)*R*V36 + 1//25*q36(t)*S*V3, q3(t)]| +|Chemical reaction network|[k6, k5, k4, k3, k2, k1, x6(t), x5(t), x4(t), x3(t), x2(t), x1(t)]|[k6, k5, k4, k3, k2, k1, x6(t), x5(t), x3(t), x1(t), x2(t), x4(t)]| +|Crauste_SI|[mu_PL, mu_PE, mu_P, mu_N, mu_M, mu_LL, mu_LE, mu_EE, delta_LM, delta_EL, rho_E//rho_P, delta_NE//rho_P, M(t), S(t), E(t), N(t), P(t)*rho_P]|[mu_PL, mu_PE, mu_P, mu_N, mu_M, mu_LL, mu_LE, mu_EE, delta_LM, delta_EL, rho_E//rho_P, delta_NE//rho_P, M(t), N(t), S(t), E(t), P(t)*rho_P]| +|Fujita|[reaction_8_k1, reaction_7_k1, reaction_6_k1, reaction_5_k2, reaction_4_k1, reaction_3_k1, reaction_2_k2, reaction_1_k1 - reaction_1_k2 - reaction_9_k1, reaction_2_k1//reaction_5_k1, a3//reaction_5_k1, a2//reaction_5_k1, a1//reaction_5_k1, pS6(t)*reaction_5_k1, pAkt_S6(t)*reaction_5_k1, S6(t)*reaction_5_k1, pAkt(t)*reaction_5_k1, Akt(t)*reaction_5_k1, pEGFR_Akt(t)*reaction_5_k1, pEGFR(t)*reaction_5_k1, EGF_EGFR(t)*reaction_5_k1*reaction_9_k1]|[reaction_8_k1, reaction_7_k1, reaction_6_k1, reaction_5_k2, reaction_4_k1, reaction_3_k1, reaction_2_k2, reaction_1_k1 - reaction_1_k2 - reaction_9_k1, reaction_2_k1//reaction_5_k1, a3//reaction_5_k1, a2//reaction_5_k1, a1//reaction_5_k1, pS6(t)*reaction_5_k1, pAkt(t)*a2 + pAkt_S6(t)*a2, pEGFR(t)*a1 + pEGFR_Akt(t)*a1, EGF_EGFR(t)*reaction_5_k1*reaction_9_k1, S6(t)*reaction_5_k1, pAkt_S6(t)*reaction_5_k1, Akt(t)*reaction_5_k1, pEGFR_Akt(t)*reaction_5_k1]| +|Goodwin oscillator|[sigma, c, b, beta + delta, beta*delta, x4(t), x1(t), (alpha*gama)//x3(t), (x2(t)*gama - x3(t)*delta)//x3(t)]|[sigma, c, b, beta + delta, beta*delta, x1(t), x4(t), (x2(t)*gama + 1//2*x3(t)*beta - 1//2*x3(t)*delta)//(alpha*gama), (alpha*gama)//x3(t)]| +|HIV|[u, h, d, b, a, lm//q, c*q^2, beta*k*q, z(t), w(t), v(t)*beta, y(t)//q, x(t)//q]|[u, h, d, b, a, lm//q, c*q^2, beta*k*q, v(t)*beta, z(t), y(t)//q, w(t), x(t)//q]| +|HIV2_io|[s, d, b, c + k1 + w1 + w2, k2*q2, c*k1 + c*w1 + c*w2 + k1*w2 + w1*w2, c*k1*w2 + c*w1*w2, k1*k2*q1 + k1*k2*q2 + k2*q2*w1, x4(t), x1(t), x3(t)*k2 - x4(t)*c, x2(t)*k1*k2 + x3(t)*k1*k2 + x3(t)*k2*w1 + x4(t)*k1*w2 + x4(t)*w1*w2]|[s, d, b, c + k1 + w1 + w2, k2*q2, c*k1 + c*w1 + c*w2 + k1*w2 + w1*w2, c*k1*w2 + c*w1*w2, k1*k2*q1 + k1*k2*q2 + k2*q2*w1, x4(t), x1(t), x3(t)*k2 - x4(t)*c, x2(t)*k1*k2 + x3(t)*k1*k2 + x3(t)*k2*w1 + x4(t)*k1*w2 + x4(t)*w1*w2]| +|HighDimNonLin|[vm, p9, p8, p7, p6, p5, p4, p3, p20, p2, p19, p18, p17, p16, p15, p14, p13, p12, p11, p10, p1, km, x20(t), x19(t), x18(t), x17(t), x16(t), x15(t), x14(t), x13(t), x12(t), x11(t), x10(t), x9(t), x8(t), x7(t), x6(t), x5(t), x4(t), x3(t), x2(t), x1(t)]|[vm, p9, p8, p7, p6, p5, p4, p3, p20, p2, p19, p18, p17, p16, p15, p14, p13, p12, p11, p10, p1, km, x20(t), x19(t), x18(t), x17(t), x16(t), x15(t), x14(t), x13(t), x12(t), x11(t), x10(t), x9(t), x8(t), x7(t), x6(t), x5(t), x4(t), x3(t), x2(t), x1(t)]| +|JAK-STAT 1|[t9, t8, t7, t6, t5, t4, t3, t20, t2, t19, t18, t16, t14, t13, t12, t10, t1, t17*t22, t15*t21, t11*t21, x10(t), x9(t), x7(t), x6(t), x5(t), x4(t), x3(t), x2(t), x1(t), x8(t)*t21]|[t9, t8, t7, t6, t5, t4, t3, t20, t2, t19, t18, t16, t14, t13, t12, t10, t1, t17*t22, t15*t21, t11*t21, x5(t), x2(t), x1(t) + x3(t) + x4(t), x9(t), (x8(t)*t17*t22)//t11, x4(t), x3(t), x7(t), x10(t), x6(t)]| +|KD1999|[Vh, V, Th, Ta, Ca0, DH//UA, E//R, (cph*roh)//UA, (cp*ro)//UA, Tj(t), T(t), Cb(t), Ca(t), Arr(t)*k0]|[Vh, V, Th, Ta, Ca0, DH//UA, E//R, (cph*roh)//UA, (cp*ro)//UA, Cb(t), Ca(t), Tj(t), T(t), Arr(t)*k0]| +|LLW1987_io|[p1 + p3, p2*p4, p1*p3, x3(t), x1(t)*x2(t), x1(t)*p4 + x2(t)*p2, x1(t)*p1*p4 + x2(t)*p2*p3]|[p1 + p3, p2*p4, p1*p3, x1(t)*p4 + x2(t)*p2, x3(t), x1(t)*x2(t), (x1(t)*p4 - x2(t)*p2)//(p1 - p3)]| +|LeukaemiaLeon2021|-|-| +|MAPK model (5 outputs bis)|-|-| +|MAPK model (5 outputs)|[gamma1110, gamma1101, gamma1100, gamma1000, gamma0100, c1011, c0111, c0011, c0010, c0001, beta11, beta10, beta01, b10, b01, b00, alpha11, alpha10, alpha01, a10, a01, a00, S11(t), S10(t), S01(t), S00(t), F(t), K(t), FS11(t), FS10(t), FS01(t), KS10(t), KS01(t), KS00(t)]|[gamma1110, gamma1101, gamma1100, gamma1000, gamma0100, c1011, c0111, c0011, c0010, c0001, beta11, beta10, beta01, b10, b01, b00, alpha11, alpha10, alpha01, a10, a01, a00, S00(t), F(t), K(t), FS11(t), S01(t), FS01(t), KS00(t), FS10(t), KS10(t), KS01(t), S11(t), S10(t)]| +|MAPK model (6 outputs)|[gamma1110, gamma1101, gamma1100, gamma1000, gamma0100, c1011, c0111, c0011, c0010, c0001, beta11, beta10, beta01, b10, b01, b00, alpha11, alpha10, alpha01, a10, a01, a00, S11(t), S10(t), S01(t), S00(t), F(t), K(t), FS11(t), FS10(t), FS01(t), KS10(t), KS01(t), KS00(t)]|[gamma1110, gamma1101, gamma1100, gamma1000, gamma0100, c1011, c0111, c0011, c0010, c0001, beta11, beta10, beta01, b10, b01, b00, alpha11, alpha10, alpha01, a10, a01, a00, S00(t), F(t), K(t), FS11(t), S01(t), FS01(t), KS00(t), FS10(t), KS10(t), KS01(t), S11(t), S10(t)]| +|Modified LV for testing|[d, a + b, a*b, x1(t), x2(t)*c]|[d, a + b, a*b, x1(t), x2(t)*c]| +|NFkB|-|-| +|PK1|[k6, k5, k4, k3 + k7, k1 + k2, k2*s3, k1*s2, k1*k3*s3, x1(t), x4(t)*s2, x3(t)*s3, x2(t)*s2]|[k6, k5, k4, k3 + k7, k1 + k2, k2*s3, k1*s2, k1*k3*s3, x1(t), x4(t)*s2, x3(t)*s3, x2(t)*s2]| +|PK2|[n, kc, ka, b2, b1, a2, a1, x3(t), x2(t), x1(t), x0(t)]|[n, kc, ka, b2, b1, a2, a1, x3(t), x1(t), x2(t), x0(t)]| +|Pharm|[n, kc, ka, b2, b1, a2, a1, x3(t), x2(t), x1(t), x0(t)]|[n, kc, ka, b2, b1, a2, a1, x3(t), x1(t), x2(t), x0(t)]| +|Phosphorylation|[k6, k5, k4, k3, k2, k1, x3(t), x1(t), x2(t), x4(t), x6(t), x5(t)]|[k6, k5, k4, k3, k2, k1, x4(t), x2(t), x1(t), x3(t), x6(t), x5(t)]| +|Pivastatin|[r3, r1, k4, k3, k2, T0*k1, T0*k, x3(t)*k1, x2(t)*k1, x1(t)*k1]|[r3, r1, k4, k3, k2, T0*k1, T0*k, x3(t)*k1, x2(t)*k1, x1(t)*k1]| +|QWWC|-|-| +|QY|[alpa, Mar, Ks, M + siga1 + siga2, phi//(Mar - siga1), M*siga1 + M*siga2 + siga1*siga2, M*siga1*siga2, (M*phi*siga2 - M*siga2)//beta_SA, (Ks*beta_SA + M*beta_SA - 1//3*Mar*beta_SA*phi + 1//3*Mar*beta_SA + 1//3*beta_SA*phi*siga2 + 2//3*beta_SA*siga1 + 2//3*beta_SA*siga2 - 1//3*beta_SI*phi*siga2 + 1//3*beta_SI*siga2)//beta_SA, P4(t), P3(t), P2(t), P1(t), P0(t), P5(t)*Mar + beta]|[alpa, Mar, Ks, M + siga1 + siga2, phi//(Mar - siga1), M*siga1 + M*siga2 + siga1*siga2, M*siga1*siga2, (M*phi*siga2 - M*siga2)//beta_SA, (Ks*beta_SA + M*beta_SA - 1//3*Mar*beta_SA*phi + 1//3*Mar*beta_SA + 1//3*beta_SA*phi*siga2 + 2//3*beta_SA*siga1 + 2//3*beta_SA*siga2 - 1//3*beta_SI*phi*siga2 + 1//3*beta_SI*siga2)//beta_SA, P3(t), P2(t), P1(t), P0(t), (P5(t)*Mar + beta)//Mar, P4(t)]| +|Ruminal lipolysis|[k4, k3, k2, x5(t), x4(t), x3(t), x2(t), x1(t)]|[k4, k3, k2, x5(t), x3(t), x1(t), x2(t), x4(t)]| +|SEAIJRC Covid model|[k, g2, g1, b, alpha, (q*r - q)//r, Ninv(t), C(t), J(t), I(t), A(t)*q, E(t)*r, S(t)*r]|[k, g2, g1, b, alpha, Ninv(t), (q*r - q)//r, C(t), J(t), I(t), A(t)*q, S(t)*r, E(t)*r]| +|SEIR 34|[mu, epsilon + gamma, epsilon*gamma, K*epsilon, beta*epsilon*r, S(t), N(t), A(t), E(t) + I(t), I(t)*gamma]|[mu, N(t), A(t), epsilon + gamma, epsilon*gamma, K*epsilon, beta*epsilon*r, I(t)*gamma, S(t), E(t) + I(t)]| +|SEIR 36 ref|[s_d, s, phi_e, phi, mu_i, mu_d, mu_0, gamma_d, gamma, beta_d, beta, F(t), Di(t), De(t), I(t), E(t), S(t), q(t), nu(t), N(t)]|[s_d, s, phi_e, phi, mu_i, mu_d, mu_0, gamma_d, gamma, beta_d, beta, q(t), nu(t), N(t), F(t), De(t), Di(t), I(t), S(t), E(t)]| +|SEIR2T|[nu, b, a, Cu(t), N(t), In(t), E(t), S(t)]|[nu, b, a, N(t), Cu(t), In(t), S(t), E(t)]| +|SEIRT|[beta, alpha + lambda, alpha*lambda, N(t), I(t), S(t)*alpha, E(t)*alpha - I(t)*lambda]|[beta, N(t), alpha + lambda, alpha*lambda, I(t), S(t)*alpha, (E(t) + I(t))//lambda]| +|SEIR_1_io|[gamma, beta//psi, gamma*psi - psi*v, psi*v - psi - v, Q(t), I(t)*psi, S(t)*psi*v, E(t)*psi*v + I(t)*psi*v]|[gamma, beta//psi, gamma*psi - psi*v, psi*v - psi - v, Q(t), I(t)*psi, S(t)*psi*v, E(t)*psi*v + I(t)*psi*v]| +|SEUIR|[z, d, (N*w)//beta, I(t), E(t)*w, S(t)*w, U(t)*w + I(t)*w]|[z, d, (N*w)//beta, I(t), U(t)*w + I(t)*w, S(t)*w, E(t)*w]| +|SIR 19|[r, q, pp, mu, beta, D(t), C(t), I(t), S(t), N(t)]|[r, q, pp, mu, beta, N(t), D(t), C(t), I(t), S(t)]| +|SIR 21|[r, q, pp, mu, beta, D(t), C(t), I(t), S(t), N(t)]|[r, q, pp, mu, beta, N(t), D(t), C(t), I(t), S(t)]| +|SIR 24|[K, gamma + mu, c*phi, I(t), S(t), A(t) - mu]|[K, gamma + mu, A(t) - mu, c*phi, I(t), S(t)]| +|SIR 6|[gamma, K//beta, N(t), I(t)*beta, S(t)*beta]|[gamma, N(t), K//beta, S(t)*beta, I(t)*beta]| +|SIRS forced|[nu, mu, g, b0, M^2, r(t), i(t), s(t), x1(t)*b1, x2(t)*M*b1]|[nu, mu, g, b0, M^2, (x2(t)*b1)//M, r(t), x1(t)*b1, i(t), s(t)]| +|SIWR original|[xi, mu, k, gam, bw, bi, a, R(t), W(t), I(t), S(t)]|[xi, mu, k, gam, bw, bi, a, W(t), R(t), I(t), S(t)]| +|SIWR with extra output|[xi, mu, k, gam, bw, bi, a, R(t), W(t), I(t), S(t)]|[xi, mu, k, gam, bw, bi, a, W(t), R(t), S(t) + I(t) + R(t), I(t)]| +|SLIQR|[s, b, Ninv, a + g, a*e*g, a*g + e*g*s - g*s, In(t), S(t)*a, Q(t)*a - Q(t)*s, L(t)*a - In(t)*g + Q(t)*s]|[s, b, Ninv, a + g, a*e*g, a*g + e*g*s - g*s, In(t), Q(t)*a - Q(t)*s, S(t)*a, L(t)*a - In(t)*g + Q(t)*s]| +|St|[T, Dd, T*a + T*d + T*dr + e + g - r - rR, (d*rR - dr*r)//(d - dr), (a*r - a*rR + d*e + d*g - dr*e - dr*g)//(d - dr), (a^2 + 2*a*d + d^2 + dr^2)//(a*dr + d*dr), (a*dr*r - a*dr*rR + d^2*g + d*dr*e - d*dr*g - dr^2*e)//(a*d - a*dr + d^2 - dr^2), S(t) + R(t), W(t)*a + W(t)*d + W(t)*dr + e + g - r - rR, S(t)*T*d - S(t)*r + R(t)*T*dr - R(t)*rR]|[T, Dd, T*a + T*d + T*dr + e + g - r - rR, (d*rR - dr*r)//(d - dr), (a*r - a*rR + d*e + d*g - dr*e - dr*g)//(d - dr), (a^2 + 2*a*d + d^2 + dr^2)//(a*dr + d*dr), (a*dr*r - a*dr*rR + d^2*g + d*dr*e - d*dr*g - dr^2*e)//(a*d - a*dr + d^2 - dr^2), S(t) + R(t), W(t)*a + W(t)*d + W(t)*dr + e + g - r - rR, (S(t)*d + R(t)*dr)//(a + d + dr)]| +|Transfection_4State|[d3, d1, b, d2//kTL, GFP(t), mRNAenz(t)*kTL, enz(t)*kTL, mRNA(t)*kTL]|[d3, d1, b, d2//kTL, GFP(t), mRNAenz(t)*kTL, enz(t)*kTL, mRNA(t)*kTL]| +|Treatment_io|[a + g + nu, b//g, d*g + nu, a*nu + g*nu, N(t), Tr(t), S(t)*g, In(t)*g - Tr(t)*nu]|[N(t), a + g + nu, b//g, d*g + nu, a*nu + g*nu, Tr(t), S(t)*g, In(t)*g - Tr(t)*nu]| +|TumorHu2019|-|-| +|TumorPillis2007|-|-| +|cLV1 (1o)|-|-| +|cLV1 (2o)|[g3, g2, g1, B31, B21, B11, A32, A31, A22, A21, A12, A11, A23//A33, A13//A33, pi2(t), pi1(t), pi3(t)*A33]|[g3, g2, g1, B31, B21, B11, A32, A31, A22, A21, A12, A11, A23//A33, A13//A33, pi2(t), pi1(t), pi3(t)*A33]| +|cholera|[xi, mu, k, gam, bw, bi, a, R(t), W(t), I(t), S(t)]|[xi, mu, k, gam, bw, bi, a, W(t), R(t), S(t) + I(t) + R(t), I(t)]| +|generalizedLoktaVolterra (1o)|[r2, r1, beta21, beta11, beta12//beta22, x1(t), x2(t)*beta22]|[r2, r1, beta21, beta11, beta12//beta22, x1(t), x2(t)*beta22]| +|generalizedLoktaVolterra (2o)|[r2, r1, beta22, beta21, beta12, beta11, x2(t), x1(t)]|[r2, r1, beta22, beta21, beta12, beta11, x2(t), x1(t)]| +|p53|[p9, p8, p7, p6, p5, p4, p3, p25, p24, p23, p21, p20, p18, p17, p16, p15, p14, p13, p12, p11, p10, p1, p22^4, x4(t), x3(t), x2(t), x1(t)]|[p9, p8, p7, p6, p5, p4, p3, p25, p24, p23, p21, p20, p18, p17, p16, p15, p14, p13, p12, p11, p10, p1, p22^4, x3(t), x2(t), x1(t), x4(t)]| *Benchmarking environment:* - - Total RAM (GiB): 15 - - Processor: Intel Xeon Processor (Icelake) - - Julia version: 1.9.3 +* Total RAM (GiB): 251 +* Processor: INTEL(R) XEON(R) GOLD 6538Y+ +* Julia version: 1.12.6 Versions of the dependencies: - - Primes : 0.5.4 - - BenchmarkTools : 1.3.2 - - IterTools : 1.8.0 - - PrecompileTools : 1.2.0 - - Symbolics : 5.5.3 - - Combinatorics : 1.0.2 - - SymbolicUtils : 1.3.0 - - DataStructures : 0.18.15 - - Groebner : 0.4.3 - - ParamPunPam : 0.0.3 - - ModelingToolkit : 8.68.0 - - AbstractAlgebra : 0.31.1 - - MacroTools : 0.5.11 - - Nemo : 0.35.3 - - SpecialFunctions : 2.3.1 +* AbstractAlgebra : 0.50.1 +* Combinatorics : 1.1.0 +* DataStructures : 0.19.6 +* Dates : 1.11.0 +* Groebner : 0.10.4 +* IterTools : 1.10.0 +* LinearAlgebra : 1.12.0 +* Logging : 1.11.0 +* MacroTools : 0.5.16 +* Nemo : 0.56.1 +* ParamPunPam : 0.5.8 +* PrecompileTools : 1.3.4 +* Primes : 0.5.7 +* Random : 1.11.0 +* RationalFunctionFields : 0.3.2 +* TimerOutputs : 0.5.29 diff --git a/benchmarking/IdentifiableFunctions/run_benchmarks.jl b/benchmarking/IdentifiableFunctions/run_benchmarks.jl index b70d03718..5c840f88e 100644 --- a/benchmarking/IdentifiableFunctions/run_benchmarks.jl +++ b/benchmarking/IdentifiableFunctions/run_benchmarks.jl @@ -1,44 +1,23 @@ #= -This script runs benchmarks and collects useful statistics. +Run identifiable-function benchmarks and compare the returned generators. -To benchmark a particular function `StructuralIdentifiability.XYZ`, execute the -following command in your favorite terminal from this directory: +From this directory, for example: +```sh +julia run_benchmarks.jl find_identifiable_functions \ + --timeout=3600 \ + --memorylimit=8192 \ + --keywords="(with_states=true, cmp=:default); (with_states=true, cmp=:cmp_lie)" ``` -$ julia run_benchmarks.jl XYZ -``` - -`XYZ` can be one of the following: -- `find_identifiable_functions`, -- `reparametrize_global`. - -Going into details, this will: -1. Load all benchmark models from benchmarks.jl and create a directory - `benchmark_results`. -2. Run the function `StructuralIdentifiability.XYZ`, collect useful runtime - statistics, and write them to `benchmark_results`. -3. Read the statistics from `benchmark_results` and produce a table with - results. - -Below in the file there is a description of some potentially useful command line -arguments. -A more involved way of running benchmarks could be the following: -``` -julia run_benchmarks.jl XYZ --target="id_total, id_are_polynomial" --timeout=3600 --keywords="(strategy=(:hybrid,),with_states=true); (strategy=(:normalforms,3),with_states=true)" -``` +Generated model files and results are written below `results/`; subprocess logs +are written below `logs/`. Both directories are ignored by Git. The summarized +comparison is written to `benchmark_result.md`. =# using ArgParse -using CpuId, Logging, Pkg, Printf -using Base.Threads -using Distributed -using Dates -using ProgressMeter - +using CpuId, Dates, Logging, Pkg, ProgressMeter using StructuralIdentifiability -using StructuralIdentifiability: _runtime_logger, ODE -using StructuralIdentifiability.ParamPunPam global_logger(Logging.ConsoleLogger(stdout, Logging.Warn)) include("../benchmarks.jl") @@ -51,7 +30,9 @@ progressbar_enabled() = include("utils.jl") -const BENCHMARK_TABLE = "benchmark_result.md" +const RESULTS_DIR = joinpath(@__DIR__, BENCHMARK_RESULTS) +const LOGS_DIR = joinpath(@__DIR__, BENCHMARK_LOGS) +const BENCHMARK_TABLE = joinpath(@__DIR__, "benchmark_result.md") function parse_commandline() s = ArgParseSettings() @@ -67,19 +48,14 @@ function parse_commandline() Semicolon-separated list of named tuples.""" arg_type = String default = "" - # default = "(strategy=(:gb, ),); (strategy=(:gb, ),with_states=true); (strategy=(:normalforms, 2),); (strategy=(:normalforms, 2),with_states=true); (strategy=(:normalforms, 3),); (strategy=(:normalforms, 3),with_states=true); (strategy=(:hybrid, ),); (strategy=(:hybrid, ),with_states=true)" - "--target" - help = """ - Statistics to be displayed in the table. - Comma-separated list of entities. - NOTE: each of these must be present in the list of - `ALL_POSSIBLE_CATEGORIES` in `utils.jl`""" - arg_type = String - default = join(ALL_CATEGORIES, ",") "--timeout" help = "Timeout, s." arg_type = Int default = 3600 + "--memorylimit" + help = "Optional memory limit in MiB for each benchmark subprocess. Use 0 to disable." + arg_type = Int + default = 0 "--workers" help = "The number of available worker processes." arg_type = Int @@ -87,25 +63,23 @@ function parse_commandline() "--skip" help = "Skip specified benchmark models." arg_type = Vector{Symbol} - default = [:NFkB,:Covid2, :Lincomp1, :Lincomp2] - # "JAK-STAT 1", - # "LeukaemiaLeon2021", - # "MAPK model (5 outputs bis)", - # "QWWC", - # "TumorHu2019", - # "TumorPillis2007", - # "cLV1 (1o)"] + default = [ + :NFkB, + :Covid2, + :Lincomp1, + :Lincomp2, + Symbol("cLV1 (1o)"), + Symbol("cLV1 (2o)"), + :TumorHu2019, + :TumorPillis2007, + ] "--models" help = """ Run specified benchmark models. - A comma-separated list of names. + A comma-separated list of benchmark IDs. Leave empty for selecting all models.""" arg_type = String default = "" - "--augment" - help = "Augment the benchmark dataset with similar models." - arg_type = Bool - default = false "--regen" help = "Re-generate the folder with benchmarks from scratch." arg_type = Bool @@ -122,71 +96,91 @@ function parse_commandline() return parse_args(s) end -function populate_benchmarks(args, kwargs) - regen = args["regen"] - dir_present = isdir((@__DIR__) * "/$BENCHMARK_RESULTS/") - dir_present && !regen && return false - @debug "Re-generating the benchmarks folder" - try - if isdir((@__DIR__) * "/$BENCHMARK_RESULTS/") - rm((@__DIR__) * "/$BENCHMARK_RESULTS/", recursive = true, force = true) - end - catch err - @info "Something went wrong when deleting the benchmarks folder" - showerror(stdout, err) - println(stdout) +function selected_benchmark_names(args) + requested_ids = filter(!isempty, strip.(split(args["models"], ","))) + names = if isempty(requested_ids) + filter(name -> isdir(joinpath(RESULTS_DIR, name)), readdir(RESULTS_DIR)) + else + [benchmarks[Symbol(id)][:name] for id in requested_ids] end - prog = Progress( + skipped_names = [ + benchmarks[id][:name] for id in args["skip"] if haskey(benchmarks, id) + ] + return sort(setdiff(names, skipped_names)) +end + +function populate_benchmarks(args) + results_present = isdir(RESULTS_DIR) + results_present && !args["regen"] && return + results_present && rm(RESULTS_DIR; recursive = true, force = true) + mkpath(RESULTS_DIR) + + @debug "Re-generating the benchmarks folder" + progress = Progress( length(benchmarks), - "Generating benchmarks", - # spinner = true, + "Generating benchmarks"; dt = 0.1, enabled = progressbar_enabled(), color = _progressbar_color, ) - for (id, bmark) in benchmarks - next!(prog) # , spinner = "⌜⌝⌟⌞") - name = bmark[:name] - system = bmark[:ode] + for (_, benchmark) in benchmarks + next!(progress) + name = benchmark[:name] + system = benchmark[:ode] @debug "Generating $name" - mkpath((@__DIR__) * "/$BENCHMARK_RESULTS/$name/") - fd = open((@__DIR__) * "/$BENCHMARK_RESULTS/$name/$name.jl", "w") - println(fd, "# $name") - println(fd, "#! format: off") - println(fd, "using StructuralIdentifiability") - println(fd, "") - ode = join(map(s -> "\t" * s, split(repr(system), "\n")[1:(end - 1)]), ",\n") - println(fd, "system = @ODEmodel(\n$ode\n)") - close(fd) + model_dir = joinpath(RESULTS_DIR, name) + mkpath(model_dir) + open(joinpath(model_dir, "$name.jl"), "w") do io + println(io, "# $name") + println(io, "#! format: off") + println(io, "using StructuralIdentifiability\n") + equations = join( + map(line -> "\t" * line, split(repr(system), "\n")[1:(end - 1)]), + ",\n", + ) + println(io, "system = @ODEmodel(\n$equations\n)") + end end - finish!(prog) - return true + finish!(progress) + return +end + +function benchmark_command(function_name, problem_name, function_kwargs, memory_limit) + project_dir = dirname(Base.active_project()) + script = joinpath(@__DIR__, "run_single_benchmark.jl") + command = `$(Base.julia_cmd()) --project=$project_dir $script $function_name $problem_name $(string(function_kwargs))` + command = Cmd(command; ignorestatus = true, detach = false, env = copy(ENV)) + memory_limit == 0 && return command + + Sys.isunix() || throw(ArgumentError("--memorylimit is supported only on Unix")) + shell_command = string( + "ulimit -v ", + memory_limit * 1024, + "; exec ", + join(Base.shell_escape.(command.exec), " "), + ) + shell = Cmd(["bash", "-lc", shell_command]) + return Cmd(shell; ignorestatus = true, detach = false, env = copy(ENV)) end function run_benchmarks(args, kwargs) - to_skip = args["skip"] timeout = args["timeout"] @assert timeout > 0 + memorylimit = args["memorylimit"] + @assert memorylimit >= 0 function_name = args["function"] @assert function_name in ("find_identifiable_functions", "reparametrize_global") nworkers = args["workers"] @assert nworkers > 0 - dirnames = first(walkdir((@__DIR__) * "/$BENCHMARK_RESULTS/"))[2] - models_from_args = filter(s -> !isempty(s), map(strip, split(args["models"], ","))) - to_run_names = if isempty(models_from_args) - dirnames - else - [benchmarks[Symbol(id)][:name] for id in models_from_args] - end - to_run_names = setdiff(to_run_names, to_skip) - to_run_indices = collect(1:length(to_run_names)) + to_run_names = selected_benchmark_names(args) + to_run_indices = eachindex(to_run_names) @info """ Benchmarking `$function_name`.""" @info """ Passing these keyword arguments to `$function_name`: - \t$(join(map(string, kwargs), "\n\t"))""" + $(join(map(string, kwargs), "\n\t"))""" @info """ Number of benchmark systems: $(length(to_run_indices)) Workers: $(nworkers) @@ -223,7 +217,6 @@ function run_benchmarks(args, kwargs) prog = Progress( length(queue), "Running benchmarks", - # spinner = true, dt = 0.3, enabled = progressbar_enabled(), color = _progressbar_color, @@ -236,21 +229,17 @@ function run_benchmarks(args, kwargs) problem_name = to_run_names[problem_id] global_id = keywords_to_global_id(function_kwargs) @debug "Running $problem_name / $global_id. Kwargs:\n$function_kwargs" - logfn = generic_filename("logs", global_id) - logs = open((@__DIR__) * "/$BENCHMARK_RESULTS/$problem_name/$logfn", "w") - # escaping ':' for Windows - function_kwargs_esc = replace(string(function_kwargs), ":" => "\\:") - cmd = Cmd( - [ - "julia", - (@__DIR__) * "/run_single_benchmark.jl", - "$function_name", - "$problem_name", - "$function_kwargs_esc", - ] + log_dir = joinpath(LOGS_DIR, problem_name) + mkpath(log_dir) + log_path = joinpath(log_dir, generic_filename("run", global_id) * ".log") + logs = open(log_path, "w") + command = benchmark_command( + function_name, + problem_name, + function_kwargs, + memorylimit, ) - cmd = Cmd(cmd, ignorestatus = true, detach = false, env = copy(ENV)) - proc = run(pipeline(cmd, stdout = logs, stderr = logs), wait = false) + proc = run(pipeline(command; stdout = logs, stderr = logs); wait = false) push!( processes, ( @@ -262,7 +251,6 @@ function run_benchmarks(args, kwargs) start_time = time_ns(), global_run_id = global_id, logfile = logs, - # errfile = errs, ), ) push!(running, processes[end]) @@ -271,7 +259,6 @@ function run_benchmarks(args, kwargs) showvalues = generate_showvalues(processes), step = 0, valuecolor = _progressbar_value_color, - # spinner = "⌜⌝⌟⌞", ) end @@ -285,7 +272,6 @@ function run_benchmarks(args, kwargs) push!(errored, proc) end close(proc.logfile) - # close(proc.errfile) kw = proc.function_kwargs start_time = proc.start_time next!( @@ -300,8 +286,9 @@ function run_benchmarks(args, kwargs) if seconds_passed(start_time) > timeout push!(to_be_removed, i) kill(proc.julia_process) + wait(proc.julia_process) close(proc.logfile) - # close(proc.errfile) + push!(errored, proc) kw = proc.function_kwargs next!( prog, @@ -321,7 +308,7 @@ function run_benchmarks(args, kwargs) finish!(prog) if !isempty(errored) - printstyled("(!) Maybe errored:\n", color = :red) + printstyled("Failed or timed out:\n", color = :red) for proc in errored println("\t$(proc.problem_name) / $(proc.global_run_id)") end @@ -330,90 +317,70 @@ function run_benchmarks(args, kwargs) return to_run_names end -function collect_timings(args, kwargs, names; content = :compare) +function compact_result_string(content) + content = strip(content) + isempty(content) && return "-" + first_bracket = findfirst(==('['), content) + if first_bracket !== nothing + content = content[first_bracket:end] + end + content = replace(content, "\n" => " ") + content = replace(content, r"\s+" => " ") + content = replace(content, "|" => "\\|") + return strip(content) +end + +function display_keyword_label(keyword) + isempty(keyword) && return "default" + if haskey(keyword, :cmp) + cmp_label = string(keyword.cmp) + if haskey(keyword, :with_states) && keyword.with_states + return "funcs ($cmp_label)" + end + return string(cmp_label) + end + return string(keywords_to_global_id(keyword)) +end + +function collect_results(args, kwargs, names) function_name = args["function"] - targets = map(Symbol, map(strip, split(args["target"], ","))) - @assert all(target -> target in ALL_CATEGORIES, targets) - @assert length(targets) > 0 @info """ Collecting benchmark results for `$function_name`. Keyword arguments of interest: - \t$(join(map(string, kwargs), "\n\t")) - - Statistics of interest: - \t$(join(map(string, targets), "\n\t")) + $(join(map(string, kwargs), "\n\t")) """ cannot_collect = [] names = sort(names) kwids = map(keywords_to_global_id, kwargs) + labels = map(display_keyword_label, kwargs) + memory_limit = args["memorylimit"] + memory_description = memory_limit == 0 ? "disabled" : "$memory_limit MiB" - # Collect timings and data from directory BENCHMARK_RESULTS. - data = Dict() + data = Dict{String, Dict{Symbol, String}}() for name in names - @debug "==== Reading $name" - data[name] = Dict() - for kwid in kwids - timings_file = nothing - ##### - data[name][kwid] = Dict() - timingsfn = timings_filename(kwid) - try - @debug "==== Opening /$BENCHMARK_RESULTS/$name/$timingsfn" - timings_file = - open((@__DIR__) * "/$BENCHMARK_RESULTS/$name/$timingsfn", "r") - catch e - @debug "Cannot collect timings for $name / $kwid" - push!(cannot_collect, (name, kwid)) - continue + model_data = Dict{Symbol, String}() + data[name] = model_data + for id in kwids + result_path = joinpath(RESULTS_DIR, name, result_filename(id)) + if isfile(result_path) + model_data[id] = compact_result_string(read(result_path, String)) + else + @debug "Cannot collect result for $name / $id" + push!(cannot_collect, (name, id)) + model_data[id] = "-" end - lines = readlines(timings_file) - if isempty(lines) - @debug "Cannot collect timings for $name / $kwid" - push!(cannot_collect, (name, kwid)) - continue - end - @assert lines[1] == name - for line in lines[2:end] - k, v = split(line, ", ") - data[name][kwid][Symbol(k)] = parse(Float64, v) - end - close(timings_file) - ##### - datafn = data_filename(kwid) - data_file = nothing - try - @debug "==== Opening /$BENCHMARK_RESULTS/$name/$datafn" - data_file = open((@__DIR__) * "/$BENCHMARK_RESULTS/$name/$datafn", "r") - catch e - @debug "Cannot collect data for $name / $kwid" - push!(cannot_collect, (name, kwid)) - continue - end - lines = readlines(data_file) - if isempty(lines) - @debug "Cannot collect data for $name / $kwid" - push!(cannot_collect, (name, kwid)) - continue - end - @assert lines[1] == name - for line in lines[2:end] - k, v = map(strip, split(line, ",")) - data[name][kwid][Symbol(k)] = v - end - close(data_file) end end if !isempty(cannot_collect) - printstyled("(!) Cannot collect timings for:\n", color = :red) + printstyled("(!) Cannot collect results for:\n", color = :red) for (name, kwid) in cannot_collect println("\t$name / $kwid") end end - # Print the table to BENCHMARK_TABLE. resulting_md = "" resulting_md *= """ ## Benchmark results @@ -423,33 +390,17 @@ function collect_timings(args, kwargs, names; content = :compare) - Benchmarked function: `$(args["function"])` - Workers: $(args["workers"]) - Timeout: $(args["timeout"]) s - - **All timings in seconds.** + - Memory limit: $memory_description """ - makecolname(kw, target) = - kw === Symbol("") ? HUMAN_READABLE_CATEGORIES[target] : - Symbol(Symbol(kw), Symbol(" / "), HUMAN_READABLE_CATEGORIES[target]) - columns = [makecolname(kwid, target) for kwid in kwids for target in targets] - resulting_md *= "|Model|" * join(map(string, columns), "|") * "|\n" - resulting_md *= "|-----|" * join(["---" for _ in columns], "|") * "|\n" + columns = ["Model"; labels] + resulting_md *= "|" * join(columns, "|") * "|\n" + resulting_md *= "|" * join(["---" for _ in columns], "|") * "|\n" for name in names - model_data = data[name] resulting_md *= "|$name|" for kwid in kwids - if !haskey(model_data, kwid) - resulting_md *= (" - " * "|")^length(columns) - continue - end - for target in targets - if !haskey(model_data[kwid], target) - resulting_md *= " - " * "|" - else - formatting_style = CATEGORY_FORMAT[target] - resulting_md *= formatting_style(model_data[kwid][target]) * "|" - end - end + resulting_md *= data[name][kwid] * "|" end resulting_md *= "\n" end @@ -461,14 +412,13 @@ function collect_timings(args, kwargs, names; content = :compare) resulting_md *= "Versions of the dependencies:\n\n" deps = Pkg.dependencies() - stid_info = deps[findfirst(x -> x.name == "StructuralIdentifiability", deps)] - for (s, uid) in stid_info.dependencies - if deps[uid].version !== nothing - resulting_md *= "* $s : $(deps[uid].version)\n" - end + dependencies = sort(collect(Pkg.project().dependencies); by = first) + for (name, uuid) in dependencies + version = deps[uuid].version + version === nothing || (resulting_md *= "* $name : $version\n") end - return open((@__DIR__) * "/$BENCHMARK_TABLE", "w") do io + return open(BENCHMARK_TABLE, "w") do io write(io, resulting_md) end end @@ -484,32 +434,27 @@ function main() @debug "Keywords for `$(args["function"])`" @debug kwargs if args["tableonly"] - dirnames = first(walkdir((@__DIR__) * "/$BENCHMARK_RESULTS/"))[2] - models_from_args = filter(s -> !isempty(s), map(strip, split(args["models"], ","))) - problems = if isempty(models_from_args) - dirnames - else - [benchmarks[Symbol(id)][:name] for id in models_from_args] - end - problems = setdiff(problems, args["skip"]) + problems = selected_benchmark_names(args) else - flag = populate_benchmarks(args, kwargs) + populate_benchmarks(args) problems = run_benchmarks(args, kwargs) printstyled( """ Benchmarking had finished in $(round((time_ns() - timestamp) / 1.0e9, digits = 2)) seconds. - Results are written to /$BENCHMARK_RESULTS + Results are written to $RESULTS_DIR """, color = :light_green, ) end - collect_timings(args, kwargs, problems) + collect_results(args, kwargs, problems) return printstyled( """ - Table with results is written to /$BENCHMARK_TABLE + Table with results is written to $BENCHMARK_TABLE """, color = :light_green, ) end -main() +if abspath(PROGRAM_FILE) == @__FILE__ + main() +end diff --git a/benchmarking/IdentifiableFunctions/run_single_benchmark.jl b/benchmarking/IdentifiableFunctions/run_single_benchmark.jl index 42dd28cf2..15cd2aa03 100644 --- a/benchmarking/IdentifiableFunctions/run_single_benchmark.jl +++ b/benchmarking/IdentifiableFunctions/run_single_benchmark.jl @@ -1,22 +1,16 @@ -using CpuId, Logging, Pkg, Printf -using Statistics - +using Logging using StructuralIdentifiability -using StructuralIdentifiability: _runtime_logger, ODE -using StructuralIdentifiability.ParamPunPam -logger = Logging.ConsoleLogger(stdout, Logging.Info) -global_logger(logger) +global_logger(Logging.ConsoleLogger(stdout, Logging.Info)) include("utils.jl") -const data = Dict() - -const NUM_RUNS = 5 +const NUM_RUNS = 1 const FUNCTION_NAME = ARGS[1] const PROBLEM_NAME = ARGS[2] -const KWARGS = parse_keywords(ARGS[3])[1] +const KWARGS = only(parse_keywords(ARGS[3])) const GLOBAL_ID = keywords_to_global_id(KWARGS) +const RESULT_DIR = joinpath(@__DIR__, BENCHMARK_RESULTS, PROBLEM_NAME) @info "" FUNCTION_NAME @info "" PROBLEM_NAME @@ -25,98 +19,79 @@ const GLOBAL_ID = keywords_to_global_id(KWARGS) flush(stdout) flush(stderr) -# Load the system -path = (@__DIR__) * "/$BENCHMARK_RESULTS/$PROBLEM_NAME/$PROBLEM_NAME.jl" -include(path) +include(joinpath(RESULT_DIR, "$PROBLEM_NAME.jl")) + +const RESOLVED_KWARGS = if haskey(KWARGS, :cmp) && KWARGS.cmp isa Symbol + merge(KWARGS, (cmp = benchmark_comparator(KWARGS.cmp, system),)) +else + KWARGS +end +const BENCHMARK_FUNCTION = getproperty(StructuralIdentifiability, Symbol(FUNCTION_NAME)) -# Compile -const FUNC_SYM = Symbol(FUNCTION_NAME) -ode = @ODEmodel(x1'(t) = a * x1(t) + x2, x2'(t) = b * x2(t) + c * d, y(t) = x1(t)) -eval(:($FUNC_SYM($system; ($KWARGS...)))) +# Compile before collecting timings. +BENCHMARK_FUNCTION(system; RESOLVED_KWARGS...) function process_system() @info "Processing $PROBLEM_NAME" @info """ Averaging over $NUM_RUNS runs. Using keyword arguments: - $(typeof(KWARGS)) - $KWARGS + $RESOLVED_KWARGS ID: $GLOBAL_ID""" - data[PROBLEM_NAME] = Dict{Any, Any}(c => 0.0 for c in ALL_CATEGORIES) + model_data = Dict{Symbol, Any}(category => 0.0 for category in ID_TIME_CATEGORIES) for _ in 1:NUM_RUNS - timing = @timed result = eval(:($FUNC_SYM($system; ($KWARGS...)))) - data[PROBLEM_NAME][:return_value] = result - @info "Result is" result - for cat in ID_TIME_CATEGORIES - if haskey(StructuralIdentifiability._runtime_logger, cat) - data[PROBLEM_NAME][cat] = StructuralIdentifiability._runtime_logger[cat] + timing = @timed BENCHMARK_FUNCTION(system; RESOLVED_KWARGS...) + model_data[:return_value] = timing.value + @info "Result is" timing.value + + for category in ID_TIME_CATEGORIES + category === :id_total && continue + if haskey(StructuralIdentifiability._runtime_logger, category) + model_data[category] += StructuralIdentifiability._runtime_logger[category] end end - for cat in ID_DATA_CATEGORIES - if haskey(StructuralIdentifiability._runtime_logger, cat) - data[PROBLEM_NAME][cat] = - deepcopy(StructuralIdentifiability._runtime_logger[cat]) + for category in ID_DATA_CATEGORIES + if haskey(StructuralIdentifiability._runtime_logger, category) + model_data[category] = + deepcopy(StructuralIdentifiability._runtime_logger[category]) end end - data[PROBLEM_NAME][:id_total] = timing.time + model_data[:id_total] += timing.time end - for cat in ID_TIME_CATEGORIES - if haskey(data[PROBLEM_NAME], cat) - data[PROBLEM_NAME][cat] = data[PROBLEM_NAME][cat] / NUM_RUNS - end + + for category in ID_TIME_CATEGORIES + model_data[category] /= NUM_RUNS end - return + return model_data end -function dump_timings() - timings = "" - timings *= "$PROBLEM_NAME\n" - for (key, model_data) in data - for c in ID_TIME_CATEGORIES - timings *= "$c, " - timings *= string(model_data[c]) * "\n" +function dump_timings(model_data) + filename = joinpath(RESULT_DIR, timings_filename(GLOBAL_ID)) + return open(filename, "w") do io + println(io, PROBLEM_NAME) + for category in ID_TIME_CATEGORIES + println(io, "$category, $(model_data[category])") end end - filename = timings_filename(GLOBAL_ID) - return open((@__DIR__) * "/$BENCHMARK_RESULTS/$PROBLEM_NAME/$filename", "w") do io - write(io, timings) - end end -function dump_results() - filename = result_filename(GLOBAL_ID) - open((@__DIR__) * "/$BENCHMARK_RESULTS/$PROBLEM_NAME/$filename", "w") do io - if haskey(data, PROBLEM_NAME) - println(io, data[PROBLEM_NAME][:return_value]) - end +function dump_results(model_data) + result_path = joinpath(RESULT_DIR, result_filename(GLOBAL_ID)) + open(result_path, "w") do io + println(io, model_data[:return_value]) end - filename = data_filename(GLOBAL_ID) - open((@__DIR__) * "/$BENCHMARK_RESULTS/$PROBLEM_NAME/$filename", "w") do io - write(io, "$PROBLEM_NAME\n") - end - for cat in ID_DATA_CATEGORIES - if !haskey(data[PROBLEM_NAME], cat) - continue - end - if cat === :something_important - # make a separate file for it - filename_cat = generic_filename(cat, GLOBAL_ID) - open((@__DIR__) * "/$BENCHMARK_RESULTS/$PROBLEM_NAME/$filename_cat", "w") do io - # print something - end - continue - end - # otherwise, print in the data file - open((@__DIR__) * "/$BENCHMARK_RESULTS/$PROBLEM_NAME/$filename", "a+") do io - write(io, "$cat, ") - write(io, string(data[PROBLEM_NAME][cat])) - write(io, "\n") + + data_path = joinpath(RESULT_DIR, data_filename(GLOBAL_ID)) + return open(data_path, "w") do io + println(io, PROBLEM_NAME) + for category in ID_DATA_CATEGORIES + haskey(model_data, category) || continue + println(io, "$category, $(model_data[category])") end end - return end -process_system() -dump_timings() -dump_results() +const MODEL_DATA = process_system() +dump_timings(MODEL_DATA) +dump_results(MODEL_DATA) diff --git a/benchmarking/IdentifiableFunctions/utils.jl b/benchmarking/IdentifiableFunctions/utils.jl index b1906b99b..a394b0b84 100644 --- a/benchmarking/IdentifiableFunctions/utils.jl +++ b/benchmarking/IdentifiableFunctions/utils.jl @@ -1,7 +1,7 @@ # shared utilities for benchmarking -using Printf const BENCHMARK_RESULTS = "results" +const BENCHMARK_LOGS = "logs" const ID_TIME_CATEGORIES = [ :id_total, @@ -18,108 +18,63 @@ const ID_TIME_CATEGORIES = [ const ID_DATA_CATEGORIES = [:id_npoints_degree, :id_npoints_interpolation, :are_id_funcs_polynomial] const ALL_CATEGORIES = union(ID_TIME_CATEGORIES, ID_DATA_CATEGORIES) -const ALL_POSSIBLE_CATEGORIES = union(ALL_CATEGORIES, Symbol[]) -const HUMAN_READABLE_CATEGORIES = Dict( - :id_io_time => "io", - :id_primality_evaluate => "io/primality-evaluate", - :id_uncertain_factorization => "io/uncertain-factor", - :id_global_time => "global id.", - :id_ideal_time => "gen. ideal", - :id_inclusion_check => "inclusion", - :id_inclusion_check_mod_p => "inclusion Zp", - :id_groebner_time => "ParamPunPam.jl", - :id_total => "Runtime", - :id_beautifulization => "beautifulization", - :id_normalforms_time => "normal forms", - :id_gbfan_time => "GB fan", - :id_ranking => "Score", - :implicit_relations => "Algebraic relations", - :dim_before => "Dim. before", - :dim_after => "Dim. after", - :are_id_funcs_polynomial => "Polynomial?", - :id_npoints_degree => "# Points, degree", - :id_npoints_interpolation => "# Points, interpolation", -) +function parse_keywords(input) + isempty(strip(input)) && return [NamedTuple()] -const CATEGORY_FORMAT = Dict() -for cat in ALL_POSSIBLE_CATEGORIES - CATEGORY_FORMAT[cat] = (val) -> if val isa Real - @sprintf("%.2f", val) - else - string(val) + return map(split(input, ";")) do keyword_set + expression = Meta.parse(strip(keyword_set)) + keywords = Core.eval(@__MODULE__, expression) + keywords == () && return NamedTuple() + keywords isa NamedTuple || + throw(ArgumentError("Expected a named tuple of keyword arguments, got $keywords")) + return keywords end end -CATEGORY_FORMAT[:are_id_funcs_polynomial] = (val) -> string(val) == "true" ? "yes" : "no" -function parse_keywords(keywords) - keywords = replace(keywords, "\\:" => ":") - keywords = replace(keywords, "\"" => "") - if isempty(keywords) - return [[]] +function benchmark_comparator(label, ode) + if label in (:default, :cmp_default) + return StructuralIdentifiability.default_cmp(ode) + elseif label === :cmp_lie + return StructuralIdentifiability.cmp_lie( + ode, + StructuralIdentifiability.RationalFunctionFields.rational_function_cmp, + ) end - if isempty(eval(Meta.parse(keywords))) - return [[]] - end - sets_of_keywords = map(strip, split(keywords, ";")) - @assert !isempty(sets_of_keywords) - kws_named_tuples = [] - for kwset in sets_of_keywords - @debug "" kwset - nt = eval(Meta.parse(kwset)) - push!(kws_named_tuples, nt) - end - return kws_named_tuples + throw(ArgumentError("Unknown comparator label: $label")) +end + +function benchmark_keyword_label(value) + value isa Symbol && return value + label = value isa Function ? typeof(value) : value + return Symbol(replace(string(label), r"[^A-Za-z0-9_]+" => "_")) end function keywords_to_global_id(keywords) - if isempty(keywords) - return Symbol() - end - id = get(keywords, :strategy, Symbol("")) - if haskey(keywords, :with_states) - if keywords.with_states - id = if id !== Symbol("") - Symbol(id, :_with_states) - else - Symbol(id, :with_states) - end - end + isempty(keywords) && return Symbol() + + id = get(keywords, :strategy, Symbol()) + if get(keywords, :with_states, false) + id = isempty(string(id)) ? :with_states : Symbol(id, :_with_states) end if haskey(keywords, :rational_interpolator) - interpolator = keywords.rational_interpolator - id = Symbol(id, :_, interpolator) + id = Symbol(id, :_, keywords.rational_interpolator) + end + if haskey(keywords, :cmp) + comparator = benchmark_keyword_label(keywords.cmp) + id = isempty(string(id)) ? comparator : Symbol(id, :_, comparator) end - if haskey(keywords, :adjoin_identifiable) - if keywords.adjoin_identifiable - id = if id !== Symbol("") - Symbol(id, :_adjoin_identifiable) - else - Symbol(id, :adjoin_identifiable) - end - end + if get(keywords, :adjoin_identifiable, false) + id = isempty(string(id)) ? :adjoin_identifiable : Symbol(id, :_adjoin_identifiable) end return id end -function timings_filename(kwid) - return generic_filename("timings", kwid) -end - -function result_filename(kwid) - return generic_filename("result", kwid) -end +timings_filename(id) = generic_filename("timings", id) +result_filename(id) = generic_filename("result", id) +data_filename(id) = generic_filename("data", id) -function data_filename(kwid) - return generic_filename("data", kwid) -end - -function generic_filename(name, kwid) - str = if kwid === Symbol("") - "$name" - else - "$(name)_$kwid" - end - str = replace(str, ":" => "") - return str +function generic_filename(name, id) + filename = isempty(string(id)) ? name : "$(name)_$id" + return replace(filename, ":" => "") end diff --git a/benchmarking/benchmark_result.md b/benchmarking/benchmark_result.md index 7eeed7d4a..86518c171 100644 --- a/benchmarking/benchmark_result.md +++ b/benchmarking/benchmark_result.md @@ -1,38 +1,28 @@ -| Model | loc_time | glob_time | ioeq_time | wrnsk_time | rank_time | check_time | total | -|:-------------------------- |:-------- |:--------- |:--------- |:---------- |:--------- |:---------- |:------- | -| Goodwin oscillator | 0.03 | 0.09 | 0.02 | 0.04 | 0.00 | 0.02 | 0.12 | -| Akt pathway | 0.64 | 4.37 | 0.18 | 0.95 | 0.00 | 3.23 | 5.01 | -| SIWR original | 0.08 | 13.08 | 2.97 | 2.86 | 0.14 | 7.10 | 13.16 | -| SEAIJRC Covid model | 0.07 | 76.34 | 27.58 | 16.23 | 0.63 | 31.90 | 76.41 | -| MAPK model (6 outputs) | 4.96 | 31.59 | 4.24 | 10.70 | 0.00 | 16.65 | 36.55 | -| Pharm | 0.05 | 100.86 | 19.30 | 39.82 | 5.64 | 36.10 | 100.90 | -| SIRS forced | 0.04 | 8.97 | 0.89 | 5.71 | 0.22 | 2.14 | 9.01 | -| SIWR with extra output | 0.09 | 0.67 | 0.44 | 0.09 | 0.00 | 0.13 | 0.76 | -| CD8 T cell differentiation | 0.48 | 0.20 | 0.02 | 0.09 | 0.00 | 0.09 | 0.69 | -| HIV | 0.11 | 0.41 | 0.04 | 0.06 | 0.00 | 0.31 | 0.52 | -| MAPK model (5 outputs bis) | 10.33 | 1665.97 | 461.87 | 214.31 | 0.05 | 989.75 | 1676.31 | -| Chemical reaction network | 0.04 | 0.14 | 0.05 | 0.05 | 0.00 | 0.04 | 0.18 | -| Modified LV for testing | 1.17 | 2.46 | 1.50 | 0.63 | 0.01 | 0.24 | 11.59 | -| MAPK model (5 outputs) | 3.91 | 98.06 | 59.62 | 9.97 | 0.00 | 28.46 | 101.97 | +|Model|loc_time|glob_time|wrnsk_time|rank_time|check_time|total| +|-----|---|---|---|---|---|---| +|SEUIR|0.61|5.91|0.17|0.00|4.78|11.21| *Benchmarking environment:* - - Total RAM (Mb): 16384.0 - - Processor: Intel(R) Core(TM) i5-8210Y CPU @ 1.60GHz - - Julia version: 1.7.2 +* Total RAM (Mb): 257647.125 +* Processor: INTEL(R) XEON(R) GOLD 6538Y+ +* Julia version: 1.12.6 Versions of the dependencies: - - Primes : 0.5.2 - - BenchmarkTools : 1.3.1 - - Singular : 0.10.1 - - IterTools : 1.4.0 - - DataStructures : 0.18.11 - - Groebner : 0.2.5 - - ModelingToolkit : 8.9.0 - - AbstractAlgebra : 0.25.3 - - MacroTools : 0.5.9 - - Nemo : 0.30.0 - - Hecke : 0.13.4 - - SpecialFunctions : 2.1.4 - - TestSetExtensions : 2.0.0 +* Random : 1.11.0 +* Primes : 0.5.7 +* TimerOutputs : 0.5.29 +* IterTools : 1.10.0 +* PrecompileTools : 1.3.4 +* LinearAlgebra : 1.12.0 +* Combinatorics : 1.1.0 +* DataStructures : 0.18.22 +* Groebner : 0.10.6 +* ParamPunPam : 0.5.8 +* RationalFunctionFields : 0.3.2 +* AbstractAlgebra : 0.46.5 +* MacroTools : 0.5.16 +* Dates : 1.11.0 +* Nemo : 0.51.1 +* Logging : 1.11.0 diff --git a/docs/src/tutorials/identifiable_functions.md b/docs/src/tutorials/identifiable_functions.md index c161af2a5..b297e6336 100644 --- a/docs/src/tutorials/identifiable_functions.md +++ b/docs/src/tutorials/identifiable_functions.md @@ -40,6 +40,16 @@ By default, `find_identifiable_functions` tries to simplify the output functions the degree of simplification. The default value is `:standard` but one could use `:strong` to try to simplify further (at the expense of heavier computation) or use `:weak` to simplify less (but compute faster). +When several generating sets describe the same field, the `cmp` keyword can be used to control which rational functions are considered simpler and end up in the final result. It is possible to provide a custom function `cmp(f, g)` that returns `true` when `f` should be preferred to `g`. For example, the following comparator prefers functions with shorter printed representations: + +```@example funcs +find_identifiable_functions( + LLW1987, + with_states = true, + cmp = (f1, f2) -> length(string(f1)) < length(string(f2)), +) +``` + As `assess_identifiability` and `assess_local_identifiability`, `find_identifiable_functions` accepts an optional parameter `loglevel` (default: `Logging.Info`) to adjust the verbosity of logging. diff --git a/ext/ModelingToolkitSIExt.jl b/ext/ModelingToolkitSIExt.jl index 85c13e15f..ee382efe0 100644 --- a/ext/ModelingToolkitSIExt.jl +++ b/ext/ModelingToolkitSIExt.jl @@ -651,11 +651,12 @@ system. ## Options -This functions takes the following optional arguments: +This function takes the following optional arguments: - `measured_quantities` - the output functions of the model. - `known_ic` - a list of functions whose initial conditions are assumed to be known, then the returned identifiable functions will be functions of parameters and initial conditions, not states (this is an experimental functionality). +- `cmp` - a comparator for rational functions. - `loglevel` - the verbosity of the logging (can be Logging.Error, Logging.Warn, Logging.Info, Logging.Debug) @@ -693,6 +694,7 @@ function StructuralIdentifiability.find_identifiable_functions( simplify = :standard, rational_interpolator = :VanDerHoevenLecerf, loglevel = Logging.Info, + cmp = nothing, ) restart_logging(loglevel = loglevel) reset_timings() @@ -706,6 +708,7 @@ function StructuralIdentifiability.find_identifiable_functions( with_states = with_states, simplify = simplify, rational_interpolator = rational_interpolator, + cmp = cmp, ) end end @@ -719,9 +722,11 @@ function _find_identifiable_functions( with_states = false, simplify = :standard, rational_interpolator = :VanDerHoevenLecerf, + cmp = nothing, ) Random.seed!(seed) ode, conversion = mtk_to_si(ode, measured_quantities) + isnothing(cmp) && (cmp = StructuralIdentifiability.default_cmp(ode)) known_ic_ = [eval_at_nemo(each, conversion) for each in known_ic] result = nothing if isempty(known_ic) @@ -732,6 +737,7 @@ function _find_identifiable_functions( seed = seed, with_states = with_states, rational_interpolator = rational_interpolator, + cmp = cmp, ) else result = StructuralIdentifiability._find_identifiable_functions_kic( @@ -741,6 +747,7 @@ function _find_identifiable_functions( prob_threshold = prob_threshold, seed = seed, rational_interpolator = rational_interpolator, + cmp = cmp, ) end result = [parent_ring_change(f, ode.poly_ring) for f in result] diff --git a/src/ODE.jl b/src/ODE.jl index dd36602e9..dcb1cb8ef 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -62,10 +62,7 @@ end #------------------------------------------------------------------------------ -function add_outputs( - ode::ODE{P}, - extra_y::Dict{String, <:RingElem}, - ) where {P <: MPolyRingElem} +function add_outputs(ode::ODE{P}, extra_y::Dict{String, <:RingElem}) where {P <: MPolyRingElem} new_var_names = vcat(collect(map(var_to_str, gens(ode.poly_ring))), collect(keys(extra_y))) new_ring, new_vars = Nemo.polynomial_ring(base_ring(ode.poly_ring), new_var_names) diff --git a/src/StructuralIdentifiability.jl b/src/StructuralIdentifiability.jl index 1a3d4d74b..6aa9fedbf 100644 --- a/src/StructuralIdentifiability.jl +++ b/src/StructuralIdentifiability.jl @@ -67,7 +67,6 @@ using RationalFunctionFields: RationalFunctionField, check_algebraicity, field_contains, - generators, simplified_generating_set using RationalFunctionFields: eval_at_dict, @@ -81,7 +80,8 @@ using RationalFunctionFields: _reduce_mod_p, is_rational_func_const, fractions_to_dennums, - check_constructive_field_membership + check_constructive_field_membership, + rational_function_cmp # defining a model export ODE, @ODEmodel, @DDSmodel @@ -93,7 +93,7 @@ export assess_local_identifiability, assess_identifiability export set_parameter_values, rename_variables # extra functionality -export find_ioequations, find_identifiable_functions +export find_ioequations, find_identifiable_functions, default_cmp # exporting to other formats export print_for_maple, print_for_DAISY, print_for_COMBOS, print_for_GenSSI diff --git a/src/discrete.jl b/src/discrete.jl index 753a9c5f7..4a73a0641 100644 --- a/src/discrete.jl +++ b/src/discrete.jl @@ -57,10 +57,7 @@ end #------------------------------------------------------------------------------ # Some functions to transform DDS's -function add_outputs( - dds::DDS{P}, - extra_y::Dict{String, <:RingElem}, - ) where {P <: MPolyRingElem} +function add_outputs(dds::DDS{P}, extra_y::Dict{String, <:RingElem}) where {P <: MPolyRingElem} return DDS{P}(add_outputs(dds.ode, extra_y)) end diff --git a/src/elimination.jl b/src/elimination.jl index 0831b011c..bc631b92a 100644 --- a/src/elimination.jl +++ b/src/elimination.jl @@ -70,8 +70,8 @@ function Bezout_matrix(f::P, g::P, var_elim::P) where {P <: MPolyRingElem} for i in 1:n for j in 1:n M[i, j] = sum( - coeffs_f[j + k + 1] * coeffs_g[i - k] - - coeffs_g[j + k + 1] * coeffs_f[i - k] for k in 0:min(i - 1, n - j) + coeffs_f[j + k + 1] * coeffs_g[i - k] - coeffs_g[j + k + 1] * coeffs_f[i - k] for + k in 0:min(i - 1, n - j) ) end end @@ -188,10 +188,7 @@ mutable struct ODEPointGenerator{P} <: PointGenerator{P} cached_points::Array{Dict{P, <:FieldElem}, 1} number_type::Type - function ODEPointGenerator{P}( - ode::ODE{P}, - big_ring::MPolyRing, - ) where {P <: MPolyRingElem} + function ODEPointGenerator{P}(ode::ODE{P}, big_ring::MPolyRing) where {P <: MPolyRingElem} prec = length(ode.x_vars) + 1 number_type = typeof(one(base_ring(big_ring))) return new(ode, big_ring, prec, Array{Dict{P, number_type}}[], number_type) @@ -215,8 +212,7 @@ function Base.iterate( initial_conditions = Dict{P, Int}(x => rand(1:sample_max) for x in gpg.ode.x_vars) input_values = Dict{P, Array{Int, 1}}( - u => [rand(1:sample_max) for _ in 1:(gpg.precision)] for - u in gpg.ode.u_vars + u => [rand(1:sample_max) for _ in 1:(gpg.precision)] for u in gpg.ode.u_vars ) @debug "Computing a power series solution" ps_solution = undef diff --git a/src/global_identifiability.jl b/src/global_identifiability.jl index 4f97226ff..f201386c0 100644 --- a/src/global_identifiability.jl +++ b/src/global_identifiability.jl @@ -91,6 +91,7 @@ The function returns a tuple containing the following: with_states::Bool = false, var_change_policy = :default, rational_interpolator = :VanDerHoevenLecerf, + cmp = default_cmp(ode), ) where {T} @info "Computing IO-equations" ioeq_time = @elapsed io_equations = @@ -148,6 +149,7 @@ The function returns a tuple containing the following: seed = 42, simplify = :standard, rational_interpolator = rational_interpolator, + cmp = cmp, ) dennums_simplified = fractions_to_dennums(no_states_simplified) # switch back the ring diff --git a/src/identifiable_functions.jl b/src/identifiable_functions.jl index 2dde62697..81eb9d6db 100644 --- a/src/identifiable_functions.jl +++ b/src/identifiable_functions.jl @@ -6,7 +6,7 @@ system. ## Options -This functions takes the following optional arguments: +This function takes the following optional arguments: - `with_states`: When `true`, also reports the identifiabile functions in the ODE states. Default is `false`. - `simplify`: The extent to which the output functions are simplified. Stronger @@ -22,30 +22,30 @@ This functions takes the following optional arguments: initial conditions, not states (this is an experimental functionality). - `prob_threshold`: A float in the range from 0 to 1, the probability of correctness. Default is `0.99`. +- `cmp`: A comparator for rational functions. For two rational functions `f1` + and `f2`, `cmp(f1, f2)` should be `true` if `f1` is simpler than `f2`. + By default, `StructuralIdentifiability.default_cmp(ode)` is used. - `seed`: The rng seed. Default value is `42`. - `loglevel` - the minimal level of log messages to display (`Logging.Info` by default) ## Example ```jldoctest -julia> using StructuralIdentifiability +using StructuralIdentifiability + +ode = @ODEmodel( + x0'(t) = -(a01 + a21) * x0(t) + a12 * x1(t), + x1'(t) = a21 * x0(t) - a12 * x1(t), + y(t) = x0(t) +) -julia> using Logging +find_identifiable_functions(ode) -julia> Logging.disable_logging(Logging.Info); +# output -julia> ode = @ODEmodel( - x0'(t) = -(a01 + a21) * x0(t) + a12 * x1(t), - x1'(t) = a21 * x0(t) - a12 * x1(t), - y(t) = x0(t) - ); - -julia> Logging.disable_logging(Logging.BelowMinLevel); - -julia> identifiable_functions = find_identifiable_functions(ode; loglevel = Logging.Error); - -julia> length(identifiable_functions) -2 +2-element Vector{AbstractAlgebra.Generic.FracFieldElem{Nemo.QQMPolyRingElem}}: + a12 + a01 + a21 + a12*a01 ``` """ @@ -58,6 +58,7 @@ function find_identifiable_functions( simplify = :standard, rational_interpolator = :VanDerHoevenLecerf, loglevel = Logging.Info, + cmp = default_cmp(ode), ) where {T <: MPolyRingElem{QQFieldElem}} restart_logging(loglevel = loglevel) reset_timings() @@ -70,6 +71,7 @@ function find_identifiable_functions( with_states = with_states, simplify = simplify, rational_interpolator = rational_interpolator, + cmp = cmp, ) else id_funcs = _find_identifiable_functions_kic( @@ -79,6 +81,7 @@ function find_identifiable_functions( seed = seed, simplify = simplify, rational_interpolator = rational_interpolator, + cmp = cmp, ) # renaming variables from `x(t)` to `x(0)` return replace_with_ic(ode, id_funcs) @@ -93,6 +96,7 @@ function _find_identifiable_functions( with_states = false, simplify = :standard, rational_interpolator = :VanDerHoevenLecerf, + cmp = default_cmp(ode), ) where {T <: MPolyRingElem{QQFieldElem}} Random.seed!(seed) @assert simplify in (:standard, :weak, :strong, :absent) @@ -113,6 +117,7 @@ function _find_identifiable_functions( prob_threshold = half_p, with_states = with_states, rational_interpolator = rational_interpolator, + cmp = cmp, ) # If simplification is needed if simplify !== :absent @@ -126,7 +131,7 @@ function _find_identifiable_functions( seed = seed, simplify = simplify, rational_interpolator = rational_interpolator, - priority_variables = [parent_ring_change(p, bring) for p in ode.parameters], + cmp = cmp, ) else id_funcs_fracs = RationalFunctionFields.dennums_to_fractions(id_funcs) diff --git a/src/known_ic.jl b/src/known_ic.jl index 83f359bce..30934ae1c 100644 --- a/src/known_ic.jl +++ b/src/known_ic.jl @@ -33,6 +33,7 @@ function _find_identifiable_functions_kic( seed = 42, simplify = :standard, rational_interpolator = :VanDerHoevenLecerf, + cmp = default_cmp(ode), ) where {T <: MPolyRingElem{Nemo.QQFieldElem}} Random.seed!(seed) @assert simplify in (:standard, :weak, :strong, :absent) @@ -45,6 +46,7 @@ function _find_identifiable_functions_kic( simplify = :absent, rational_interpolator = rational_interpolator, seed = seed, + cmp = cmp, ) id_funcs = vcat(id_funcs_general, [f // one(parent(ode)) for f in known_ic]) @@ -56,6 +58,7 @@ function _find_identifiable_functions_kic( seed = seed, simplify = simplify, rational_interpolator = rational_interpolator, + cmp = cmp, ) end @@ -91,7 +94,12 @@ function _assess_identifiability_kic( funcs_to_check = vcat(ode.x_vars, ode.parameters) end half_p = 0.5 + prob_threshold / 2 - id_funcs = _find_identifiable_functions_kic(ode, known_ic, prob_threshold = half_p, simplify = :absent) + id_funcs = _find_identifiable_functions_kic( + ode, + known_ic, + prob_threshold = half_p, + simplify = :absent, + ) #funcs_to_check = replace_with_ic(ode, funcs_to_check) result = OrderedDict(f => :globally for f in funcs_to_check) diff --git a/src/output_saturation.jl b/src/output_saturation.jl index 1ccb739ba..0d0c277fb 100644 --- a/src/output_saturation.jl +++ b/src/output_saturation.jl @@ -17,18 +17,29 @@ function saturate_outputs( orders::Dict{P, Int}, max_deg = 5, ) where {P <: MPolyRingElem} - lie_ders = lie_derivatives_up_to(ode, Dict(ode.y_equations[y] => ord for (y, ord) in orders)) + lie_ders = + lie_derivatives_up_to(ode, Dict(ode.y_equations[y] => ord for (y, ord) in orders)) @info "Degrees: $([(total_degree(numerator(f)), total_degree(denominator(f))) for f in lie_ders])" - current_y = RationalFunctionField(vcat(collect(values(ode.y_equations)), ode.u_vars, ode.parameters)) + current_y = RationalFunctionField( + vcat(collect(values(ode.y_equations)), ode.u_vars, ode.parameters), + ) lie_ders = filter(f -> total_degree_frac(f) <= max_deg, lie_ders) sort!(lie_ders, lt = RationalFunctionFields.rational_function_cmp) for f in lie_ders if !first(RationalFunctionFields.check_algebraicity_modp(current_y, [f])) - current_y = RationalFunctionField(vcat(generators(current_y), [f])) + current_y = RationalFunctionField( + vcat(RationalFunctionFields.generators(current_y), [f]), + ) end end - new_outputs = generators(current_y)[(length(ode.y_vars) + length(ode.u_vars) + length(ode.parameters) + 1):end] + new_outputs = RationalFunctionFields.generators(current_y)[ + ( + length(ode.y_vars) + length( + ode.u_vars, + ) + length(ode.parameters) + 1 + ):end, + ] idx = 1 old_y_names = map(var_to_str, ode.y_vars) @@ -59,9 +70,11 @@ Inputs: Output: dictionary from outputs to orders """ function propose_orders(ode) - x_with_u = filter(x -> length(intersect(vars(ode.x_equations[x]), ode.u_vars)) > 0, ode.x_vars) + x_with_u = + filter(x -> length(intersect(vars(ode.x_equations[x]), ode.u_vars)) > 0, ode.x_vars) - graph = Dict(x => Set(intersect(vars(ode.x_equations[x]), ode.x_vars)) for x in ode.x_vars) + graph = + Dict(x => Set(intersect(vars(ode.x_equations[x]), ode.x_vars)) for x in ode.x_vars) result = Dict(y => 0 for y in ode.y_vars) for y in ode.y_vars diff --git a/src/pb_representation.jl b/src/pb_representation.jl index 96c2efa8d..756d11a10 100644 --- a/src/pb_representation.jl +++ b/src/pb_representation.jl @@ -109,7 +109,9 @@ function common_ring(poly::MPolyRingElem, pbr::PBRepresentation) "$(u)_$h" for h in 0:max( max_ords[u], - max_offset + max([difforder(p, u) for p in values(pbr.projections)]...), + max_offset + max( + [difforder(p, u) for p in values(pbr.projections)]..., + ), ) ], ) diff --git a/src/primality_check.jl b/src/primality_check.jl index 6162f621f..f085f2ff7 100644 --- a/src/primality_check.jl +++ b/src/primality_check.jl @@ -26,7 +26,8 @@ function check_primality_zerodim(J::Array{QQMPolyRingElem, 1}) @debug "Prime is $p" F = Nemo.GF(p) S_F = Nemo.matrix_space(F, dim, dim) - generic_multiplication_modp = S_F([generic_multiplication[i, j] for i in 1:dim for j in 1:dim]) + generic_multiplication_modp = + S_F([generic_multiplication[i, j] for i in 1:dim for j in 1:dim]) R, t = Nemo.polynomial_ring(F, "t") chpoly = Nemo.charpoly(R, generic_multiplication_modp) @debug "Charpoly computed" diff --git a/src/states.jl b/src/states.jl index ae48e773d..2dd7d7775 100644 --- a/src/states.jl +++ b/src/states.jl @@ -69,7 +69,10 @@ end #------------------------------------------------------------------------------ -function lie_derivatives_up_to(ode::ODE, orders::Dict{P, Int}) where {P <: AbstractAlgebra.RingElem} +function lie_derivatives_up_to( + ode::ODE, + orders::Dict{P, Int}, + ) where {P <: AbstractAlgebra.RingElem} result = Array{Generic.FracFieldElem, 1}() for (f, ord) in orders curr = extract_coefficients_ratfunc(f, ode.u_vars) @@ -117,5 +120,8 @@ identifiable functions of parameters only end end - return lie_derivatives_up_to(ode, Dict(ode.y_equations[y] => ord for (y, ord) in y_to_ord)) + return lie_derivatives_up_to( + ode, + Dict(ode.y_equations[y] => ord for (y, ord) in y_to_ord), + ) end diff --git a/src/util.jl b/src/util.jl index 81c85713b..46cdf095d 100644 --- a/src/util.jl +++ b/src/util.jl @@ -302,3 +302,71 @@ function replace_with_ic(ode, funcs) end # ----------------------------------------------------------------------------- + +""" + default_cmp(ode) + +Returns a comparator for rational functions. +Used as a default measure of simplicity in StructuralIdentifiability.jl. + +## Example + +```jldoctest; setup = :(using Logging; Logging.disable_logging(Logging.Info);) +julia> using StructuralIdentifiability + +julia> ode = @ODEmodel( + x0'(t) = -(a01 + a21) * x0(t) + a12 * x1(t), + x1'(t) = a21 * x0(t) - a12 * x1(t), + y(t) = x0(t) + ); + +julia> f1, f2 = (a01 + a21) // a12, x0 // a12; + +julia> cmp = default_cmp(ode); + +julia> cmp(f1, f2) # returns true, so f1 is considered simpler than f2 +true +``` + +""" +default_cmp(ode) = cmp_prefer_params(ode, rational_function_cmp) + +""" + cmp_prefer_params(ode, cmp) + +Wrap a comparator of rational functions `cmp` so that functions containing +only parameters of the `ode` are preferred over functions containing states. +""" +function cmp_prefer_params(ode, cmp) + params = ode.parameters + return function (f, g) + isempty(params) && return cmp(f, g) + params = map(p -> parent_ring_change(p, parent(numerator(f))), params) + prefer_f = issubset(vars(f), params) + prefer_g = issubset(vars(g), params) + prefer_f != prefer_g && return prefer_f + return cmp(f, g) + end +end + +""" + cmp_lie(ode, cmp) + +Wrap a comparator of rational functions `cmp` so Lie derivatives are compared +when they are nonzero. +""" +function cmp_lie(ode, cmp) + return function (f, g) + f = parent_ring_change(f, parent(ode)) + g = parent_ring_change(g, parent(ode)) + df = lie_derivative(f, ode) + dg = lie_derivative(g, ode) + return if iszero(df) && iszero(dg) + cmp(f, g) + else + cmp(df, dg) + end + end +end + +# ----------------------------------------------------------------------------- diff --git a/src/wronskian.jl b/src/wronskian.jl index 244077783..b2e2a86e0 100644 --- a/src/wronskian.jl +++ b/src/wronskian.jl @@ -24,9 +24,7 @@ function monomial_compress(io_equation, params::Array{<:MPolyRingElem, 1}) params_xs = isempty(params) ? empty(params) : gens(parent(first(params))) # Pre-compute param string names for faster lookup param_names = Set(var_to_str(p, xs = params_xs) for p in params) - other_vars = [ - v for v in gens(parent(io_equation)) if !(var_to_str(v) in param_names) - ] + other_vars = [v for v in gens(parent(io_equation)) if !(var_to_str(v) in param_names)] coeffdict = extract_coefficients(io_equation, other_vars) expvect = collect(keys(coeffdict)) coeffs = collect(values(coeffdict)) diff --git a/test/bodies/identifiable_functions.jl b/test/bodies/identifiable_functions.jl index 24cce99a2..237e4af7c 100644 --- a/test/bodies/identifiable_functions.jl +++ b/test/bodies/identifiable_functions.jl @@ -1113,10 +1113,10 @@ ident_funcs_states = [ k02 + k12, b1 * c, Km * c, + Vm * a1 * c + Vm * a2 * c, + (a2 * k21) // (a1 + a2), x1 * c, x2 * c * k21, - (Vm * a2) // x2, - Vm * a1 * c + Vm * a2 * c, ] push!(test_cases, (ode = ode, ident_funcs = ident_funcs, with_states = false)) push!(test_cases, (ode = ode, ident_funcs = ident_funcs_states, with_states = true)) @@ -1197,5 +1197,56 @@ push!(test_cases, (ode = ode, ident_funcs = ident_funcs_states, with_states = tr end end +@testset "Identifiable functions with custom comparator" begin + ode = @ODEmodel( + S'(t) = -beta_W * S(t) * W(t) - beta_I * S(t) * I(t), + I'(t) = beta_W * S(t) * W(t) + beta_I * S(t) * I(t) - gamma * I(t), + W'(t) = alpha * I(t) - zeta * W(t), + y(t) = W(t) + ) + + cmp1 = (f, g) -> begin + return length(string(f)) < length(string(g)) + end + + cmp2 = StructuralIdentifiability.cmp_lie( + ode, RationalFunctionFields.rational_function_cmp + ) + + cmp3 = StructuralIdentifiability.cmp_prefer_params(ode, cmp1) + + o = one(parent(ode)) + @test cmp3(beta_W // o, S // o) + @test !cmp3(S // o, beta_W // o) + + with_cmp1 = find_identifiable_functions(ode, with_states = true, cmp = cmp1) + + with_cmp2 = find_identifiable_functions(ode, with_states = true, cmp = cmp2) + + with_cmp3 = find_identifiable_functions(ode, with_states = true, cmp = cmp3) + + default_cmp = find_identifiable_functions(ode, with_states = true) + + @test fields_equal( + RationalFunctionField(with_cmp1), + RationalFunctionField(default_cmp), + 0.99, + ) + @test fields_equal( + RationalFunctionField(with_cmp2), + RationalFunctionField(default_cmp), + 0.99, + ) + @test fields_equal( + RationalFunctionField(with_cmp3), + RationalFunctionField(default_cmp), + 0.99, + ) + + @test issorted(with_cmp1, lt = cmp1) + @test issorted(with_cmp2, lt = cmp2) + @test issorted(with_cmp3, lt = cmp3) +end + # restoring the default StructuralIdentifiability.RationalFunctionFields.SAT_FACTORIZATION_DEFAULT = :none