Abstract the saturation variable index

[sumiya11]

Needs some benchmarks

[codecov]

Codecov Report

Merging #222 (fbe77de) into master (d84e7fb) will decrease coverage by 0.07%.
The diff coverage is 90.32%.

@@            Coverage Diff             @@
##           master     #222      +/-   ##
==========================================
- Coverage   90.93%   90.86%   -0.07%     
==========================================
  Files          24       24              
  Lines        3133     3143      +10     
==========================================
+ Hits         2849     2856       +7     
- Misses        284      287       +3     

Files	Coverage Δ
...rc/RationalFunctionFields/RationalFunctionField.jl	89.72% <100.00%> (ø)
src/RationalFunctionFields/normalforms.jl	97.64% <100.00%> (-0.01%)	⬇️
src/RationalFunctionFields/IdealMQS.jl	79.01% <92.30%> (-0.50%)	⬇️
src/util.jl	87.41% <84.61%> (-0.27%)	⬇️

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

[pogudingleb]

Below are the benchmarks (first and last refer to the position of the saturation variable). Most of the differences are not very significant and go different ways except for QY where first is clearly better.

Model	First	Last	First (states)	Last (states)
Akt pathway	4.23	3.86	12.98	14.61
Bilirubin2_io	2.87	3.13	16.99	18.61
Biohydrogenation_io	0.73	0.77	1.67	2.12
Bruno2016	0.63	0.67	1.23	1.55
CD8 T cell differentiation	1.33	1.19	2.64	3.63
CGV1990	4.14	3.78	7.75	8.55
Chemical reaction network	0.68	0.64	1.2	1.3
Crauste_SI	1.43	1.17	2.94	3.16
Fujita	3.82	3.43	11.2	11.53
Goodwin oscillator	0.89	0.75	2.51	2.01
HIV	1.16	1.01	3.43	2.9
HIV2_io	3.33	3.13	17.18	16.61
HighDimNonLin	52.91	56.89	64.69	59.2
JAK-STAT 1	16.26	16.26	23.99	-
KD1999	1.44	1.52	3.14	3.65
LLW1987_io	0.42	0.42	1.11	1.21
LeukaemiaLeon2021	-	-	-	-
MAPK model (5 outputs bis)	-	-	-	-
MAPK model (5 outputs)	43.23	44.6	47.81	51
MAPK model (6 outputs)	8.6	8.78	13.02	13.15
Modified LV for testing	0.32	0.3	0.62	0.72
PK1	1.38	1.34	3.53	4.31
PK2	288.78	290.42	293.44	297.41
Pharm	288.56	290.07	295.44	292.02
Phosphorylation	0.79	0.75	1.4	1.25
Pivastatin	13.32	12.52	13.46	12.71
QWWC	-	-	-	-
QY	5.47	10.42	15.27	84.42
Ruminal lipolysis	0.33	0.29	0.82	0.9
SEAIJRC Covid model	149.59	152.69	159.85	160.05
SEIR 34	0.72	0.91	2.57	2.3
SEIR 36 ref	1.82	2.15	4.05	3.86
SEIR2T	0.4	0.35	0.88	0.96
SEIRT	0.35	0.35	1.53	1.56
SEIR_1_io	0.73	0.62	2.33	2.33
SEUIR	0.7	0.71	1.69	1.51
SIR 19	0.63	0.65	1.17	1.24
SIR 21	0.67	0.61	1.19	1.17
SIR 24	0.59	0.62	1.23	1.23
SIR 6	0.1	0.07	0.82	1.22
SIRS forced	14.68	14.65	17.08	16.29
SIWR original	26.9	27.56	23.56	27.01
SIWR with extra output	1.12	1.09	1.84	1.67
SLIQR	1.64	1.58	4.32	4.63
St	50.54	52.02	142.24	145.04
Transfection_4State	0.48	0.48	1.55	1.19
Treatment_io	0.87	0.93	2.62	2.93
TumorHu2019	-	-	-	-
TumorPillis2007	-	-	-	-
cLV1 (1o)	-	-	-	-
cLV1 (2o)	2.13	1.78	4.76	4.01
cholera	1.01	1.09	1.86	1.91
generalizedLoktaVolterra (1o)	0.51	0.6	0.96	0.98
generalizedLoktaVolterra (2o)	0.62	0.54	0.97	0.84
p53	3.25	3.15	6.74	6.94

[sumiya11]

Thanks! It looks like JAK-STAT 1 also finishes with First

[sumiya11]

I think for the sake of completeness it would be interesting to run benchmarks where the saturation polynomial is factored

[sumiya11]

But that would require some tweaking in the MQS ideal. We can probably merge this with First for now

[pogudingleb]

For the record, the original model which stimulated this discussion we the following:

ode = @ODEmodel(
    x0'(t) = - D0 * x0(t),
    x1'(t) = x0(t) - D1 * x1(t),
    x2'(t) = m1 * x1(t) - D2 * x2(t),
    x3'(t) = m2 * x2(t) - D3 * x3(t), 
    x4'(t) = m3 * x3(t) - D4*x4(t),
    x02'(t) = - D02 * x02(t),
    x12'(t) = x02(t) - g12*x12(t) - (m12 + f1/(1 + z1(t))) * ( ((z2(t))^20)/(tau+(z2(t))^20) ) * x12(t),
    x22'(t) = (m12 + f1/(1 + z1(t))) * ( ((z2(t))^20)/(tau+(z2(t))^20) ) * x12(t) - D22 * x22(t),
    x32'(t) = m22 * x22(t) - D32 * x32(t),
    x42'(t) = m32 * x32(t) - D42 * x42(t),
    z1'(t) = n*z1(t)*(m22*x22(t) + (m32- D32)*x32(t)- D42*x42(t) + m2*x2(t) + (m3-D3)*x3(t) - D4*x4(t))/(x32(t)+x42(t)+x3(t)+x4(t)),
    z2'(t) = 1,
    X1'(t) = x1(t) + x12(t),
    y1(t) = x1(t) + x12(t),
    y2(t) = x2(t) + x22(t),
    y3(t) = x3(t) + x32(t),
    y4(t) = x4(t) + x42(t),
    Y1(t) = X1(t),
    Y2(t) = C2 - (X1(t)*g12/D22 + (x0(t)/D02 +x12(t) + x22(t))/D22 + x2(t)/D2 + (x1(t) + x0(t)/D0)*(D2*g12+m1*D22)/(D1*D2*D22) ),
    Y3(t) = C3 - (X1(t)*g12*m22/(D22*D32) + m22*(x0(t)/D02 + x12(t) + x22(t) + D22*x32(t)/m22)/(D22*D32) + m2*(x2(t) + x3(t)*D2/m2)/(D2*D3) + (x1(t) + x0(t)/D0)*(g12*m22*D2*D3+m1*m2*D22*D32)/(D1*D2*D3*D22*D32) ),
    Y4(t) = C4 - (m22*m32*(X1(t)*g12 + x0(t)/D02 + x12(t) + x22(t) + D22*x32(t)/m22 + D22*D32*x42(t)/(m22*m32))/(D22*D32*D42) + m3*(m2*x2(t)/D2 + x3(t) + D3*x4(t)/m3)/(D3*D4) + (x1(t) + x0(t)/D0)*(g12*m22*m32*D2*D3*D4+m1*m2*m3*D22*D32*D42)/(D1*D2*D3*D4*D22*D32*D42) ),
    Y5(t) = z2(t)
)

from this paper.

Here is also a mathematical rationale why putting the saturation variable first should be beneficial. Consider a subfield generated by

(a^3 - b^3) / (a^2 - b^2), (a^4 - b^4) / (a^2 - b^2), (a^5 - b^5) / (a^2 - b^2)

These functions are symmetric and in fact generate the field of all symmetric functions. The corresponding MQS ideal in Q(a, b)[A, B, t] is generated by

(A^4-B^4)*(a^2-b^2)-(A^2-B^2)*(a^4-b^4), (A^2-B^2)*(a^5-b^5)-(A^5-B^5)*(a^2-b^2), (A^3-B^3)*(a^2-b^2)-(A^2-B^2)*(a^3-b^3), t*(A^2-B^2)-1

Now if we compute the Groebner basis in degrevlex with A, B, t, we get:

2*B+(a^3-a^2*b-a*b^2+b^3)*t-a-b, 2*A+(-a^3+a^2*b+a*b^2-b^3)*t-a-b, (a^4-2*a^2*b^2+b^4)*t^2-1

Here we see a complicated quadratic equation on t, and then A and B are expressed via t by complicated formulas. One would have to reconstruct at least to degree 3 to get the right generators. On the other hand, let us compute Groebner basis with t, A, B. It would be

A+B-a-b, 2*B+(a^3-a^2*b-a*b^2+b^3)*t-a-b, B^2+(-a-b)*B+a*b

Here the role of "pivotal" element is taken by B and it yields much simpler expressions. So the idea is that the true parameters have physical meaning, so they will likely yield simpler expressions compared to t which is, in a way, an artefact of the original (often horrible) generators.