remove option for sparsity in Jacobian (requires a special package)

Author: TorkelE

docs/src/steady_state_functionality/steady_state_stability_computation.md

@@ -53,11 +53,5 @@ stability_2 = steady_state_stability(steady_states_2, sa_loop, ps_2; ss_jac=ss_j
 nothing # hide
 ```
 
-It is possible to designate that a [sparse Jacobian](@ref ref) should be used using the `sparse = true` option (either to `steady_state_jac` or directly to `steady_state_stability`):
-```@example stability_1
-ss_jac = steady_state_jac(sa_loop; sparse = true)
-nothing # hide
-```
-
 !!! warn
     For systems with [conservation laws](@ref homotopy_continuation_conservation_laws), `steady_state_jac` must be supplied a `u0` vector (indicating species concentrations for conservation law computation). This is required to eliminate the conserved quantities, preventing a singular Jacobian. These are supplied using the `u0` optional argument.

src/steady_state_stability.jl

@@ -2,7 +2,7 @@
 
 """
     stability(u::Vector{T}, rs::ReactionSystem, p; tol = 10*sqrt(eps())
-                sparse = false, ss_jac = steady_state_jac(u, rs, p; sparse=sparse))
+                ss_jac = steady_state_jac(u, rs, p))
 
 Compute the stability of a steady state (Returned as a `Bool`, with `true` indicating stability).
 
@@ -16,8 +16,7 @@ stable if the corresponding Jacobian's maximum eigenvalue real part is < 0. Howe
 eigenvalue is in the range `-tol< eig < tol`, and error is throw, as we do not deem that stability
 can be ensured with enough certainty. The choice `tol = 10*sqrt(eps())` has *not* been subject
 to much analysis.
-- `sparse = false`: If we wish to create a sparse Jacobian for the stability computation.
-- `ss_jac = steady_state_jac(u, rs; sparse = sparse)`: It is possible to pre-compute the 
+- `ss_jac = steady_state_jac(u, rs)`: It is possible to pre-compute the 
 Jacobian used for stability computation using `steady_state_jac`. If stability is computed 
 for many states, precomputing the Jacobian may speed up evaluation.
 
@@ -47,7 +46,7 @@ however, have not been subject to further analysis (and can be changed through t
 ```
 """
 function steady_state_stability(u::Vector, rs::ReactionSystem, ps; tol = 10*sqrt(eps(ss_val_type(u))),
-                                sparse = false, ss_jac = steady_state_jac(rs; u0 = u, sparse = sparse))
+                                ss_jac = steady_state_jac(rs; u0 = u))
     # Warning checks.
     if !is_autonomous(rs) 
         error("Attempting to compute stability for a non-autonomous system (e.g. where some rate depend on $(rs.iv)). This is not possible.")
@@ -62,6 +61,11 @@ function steady_state_stability(u::Vector, rs::ReactionSystem, ps; tol = 10*sqrt
 
     # Computes stability (by checking that the real part of all eigenvalues are negative).
     # Here, `ss_jac` is a `ODEProblem` with dummy values for `u0` and `p`.
+
+    if isdefined(Main, :Infiltrator)
+        Main.infiltrate(@__MODULE__, Base.@locals, @__FILE__, @__LINE__)
+      end
+
     J = zeros(length(u), length(u))
     ss_jac = remake(ss_jac; u0 = u, p = ps)
     ss_jac.f.jac(J, ss_jac.u0, ss_jac.p, Inf)
@@ -79,15 +83,14 @@ ss_val_type(u::Vector{Pair{S,T}}) where {S,T} = T
 ss_val_type(u::Dict{S,T}) where {S,T} = T
 
 """
-    steady_state_jac(rs::ReactionSystem; u0 = [], sparse = false)
+    steady_state_jac(rs::ReactionSystem; u0 = [])
 
 Creates the Jacobian function which can be used as input to `steady_state_stability`. Useful when 
 a large number of stability computation has to be carried out in a performant manner.
 
 Arguments:
     - `rs`: The reaction system model for which we want to compute stability.
     - `u0 = []`: For systems with conservation laws, a `u` is required to compute the conserved quantities.
-    - `sparse = false`: If we wish to create a sparse Jacobian for the stability computation.
  
 Example:
 ```julia
@@ -110,7 +113,7 @@ Notes:
 ```
 """
 function steady_state_jac(rs::ReactionSystem; u0 = [sp => 0.0 for sp in unknowns(rs)], 
-                          sparse = false, combinatoric_ratelaws = get_combinatoric_ratelaws(rs))
+                          combinatoric_ratelaws = get_combinatoric_ratelaws(rs))
     # If u0 is a vector of values, must be converted to something MTK understands.
 
     # Converts u0 to values MTK understands, and checks that potential conservation laws are accounted for.
@@ -119,7 +122,7 @@ function steady_state_jac(rs::ReactionSystem; u0 = [sp => 0.0 for sp in unknowns
 
     # Creates an `ODEProblem` with a Jacobian. Dummy values for `u0` and `ps` must be provided.
     ps = [p => 0.0 for p in parameters(rs)]
-    return ODEProblem(rs, u0, 0, ps; jac = true, remove_conserved = true, sparse = sparse, 
+    return ODEProblem(rs, u0, 0, ps; jac = true, remove_conserved = true, 
                       combinatoric_ratelaws = combinatoric_ratelaws)
 end
 

test/miscellaneous_tests/stability_computation.jl

@@ -12,7 +12,6 @@ rng = StableRNG(12345)
 
 # Tests that stability is correctly assessed (using simulation) in multi stable system.
 # Tests that `steady_state_jac` function works.
-# Tests with and without sparsity.
 # Tests using symbolic input.
 let
     # System which may have between 1 and 7 fixed points.
@@ -22,31 +21,28 @@ let
         d, (X,Y) --> 0
     end
     ss_jac = steady_state_jac(rn)
-    ss_jac_sparse = steady_state_jac(rn; sparse = true)
 
     # Repeats several times, most steady state stability cases should be encountered several times.
-    for repeat = 1:50
+    for repeat = 1:20
         # Generates random parameter values (which can generate all steady states cases).
         ps = (:v => 1.0 + 3*rand(rng), :K => 0.5 + 2*rand(rng), :n => rand(rng,[1,2,3,4]), 
               :d => 0.5 + rand(rng))
 
         # Computes stability using various jacobian options.
         sss = hc_steady_states(rn, ps)
         stabs_1 = [steady_state_stability(ss, rn, ps) for ss in sss]
-        stabs_2 = [steady_state_stability(ss, rn, ps; sparse = true) for ss in sss]
-        stabs_3 = [steady_state_stability(ss, rn, ps; ss_jac = ss_jac) for ss in sss]
-        stabs_4 = [steady_state_stability(ss, rn, ps; ss_jac = ss_jac_sparse) for ss in sss]
+        stabs_2 = [steady_state_stability(ss, rn, ps; ss_jac = ss_jac) for ss in sss]
 
         # Confirms stability using simulations.
         for (idx, ss) in enumerate(sss)
             ssprob = SteadyStateProblem(rn, [1.001, 0.999] .* ss, ps)
             sol = solve(ssprob, DynamicSS(Vern7()); abstol = 1e-8, reltol = 1e-8)
-            stabs_5 = isapprox(ss, sol.u; atol = 1e-6)
-            @test stabs_1[idx] == stabs_2[idx] == stabs_3[idx] == stabs_4[idx] == stabs_5
+            stabs_3 = isapprox(ss, sol.u; atol = 1e-6)
+            @test stabs_1[idx] == stabs_2[idx] == stabs_3
 
             # Checks stability when steady state is given on a pair form ([X => x_val, Y => y_val]).
-            stabs_6 = steady_state_stability(Pair.(unknowns(rn), ss), rn, ps)
-            @test stabs_5 == stabs_6
+            stabs_4 = steady_state_stability(Pair.(unknowns(rn), ss), rn, ps)
+            @test stabs_3 == stabs_4
         end
     end
 end