Merge pull request #19 from palmtree2013/main

adapt the code
voduchuy · Oct 22, 2022 · 49b895a · 49b895a
2 parents d2e4220 + d2a32f1
commit 49b895a
Show file tree

Hide file tree

Showing 5 changed files with 16 additions and 10 deletions.
diff --git a/Project.toml b/Project.toml
@@ -19,7 +19,7 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
-Catalyst = "10"
+Catalyst = "10, 11, 12"
 DataStructures = "0.18, 0.19"
 DifferentialEquations = "7"
 DocStringExtensions = "0.8"

diff --git a/src/cmemodel/cmemodel.jl b/src/cmemodel/cmemodel.jl
@@ -174,7 +174,7 @@ end
 
 Construct a CME model with sensitivity information using a net stoichiometry matrix `stoich_matrix` and propensity functions `propensities` and parameters `parameters`. The partial derivatives of the propensity functions with respect to parameters are automatically derived using automatic differentiation.
 """
-CmeModelWithSensitivity(stoich_matrix, propensities::Vector{<:Propensity}, parameters::AbstractVector) where {PT <: Propensity} = CmeModelWithSensitivity(CmeModel(stoich_matrix, propensities, parameters))
+CmeModelWithSensitivity(stoich_matrix, propensities::Vector{<:Propensity}, parameters::AbstractVector) = CmeModelWithSensitivity(CmeModel(stoich_matrix, propensities, parameters))
 
 # Pretty printing for CmeModelWithSensitivity
 using Printf 

diff --git a/src/cmemodel/senstools/sparsity_pattern.jl b/src/cmemodel/senstools/sparsity_pattern.jl
@@ -2,7 +2,7 @@ export propensitygrad_sparsity_pattern
 
 function propensitygrad_sparsity_pattern(speciescount, parametercount, propensities, parameters)    
     propensitycount = length(propensities)
-    ModelingToolkit.@variables t, x[1:speciescount](t), alphas[1:propensitycount](t)
+    ModelingToolkit.@variables t, (x(t))[1:speciescount], (alphas(t))[1:propensitycount]
     ModelingToolkit.@parameters p[1:parametercount]
     eqs = Vector{ModelingToolkit.Equation}([])
     for r in 1:propensitycount 

diff --git a/src/transientcme/sparse/fspsolve.jl b/src/transientcme/sparse/fspsolve.jl
@@ -11,7 +11,7 @@ function solve(model::CmeModel,
     initial_distribution::FspVectorSparse{NS,IntT,RealT},
     tspan::Union{Vector,Tuple},
     ode_method::Union{Nothing,AbstractODEAlgorithm}; 
-    saveat = [], fsptol::AbstractFloat = 1.0E-6, odeatol::AbstractFloat = 1.0E-10, odertol::AbstractFloat = 1.0E-4) where {NS,IntT<:Integer,RealT<:AbstractFloat}
+    saveat = [], fsptol::AbstractFloat = 1.0E-6, odeatol::AbstractFloat = 1.0E-6, odertol::AbstractFloat = 1.0E-4) where {NS,IntT<:Integer,RealT<:AbstractFloat}
 
     statespace = StateSpaceSparse(model.stoich_matrix, initial_distribution.states)
     sink_count = get_sink_count(statespace)
@@ -23,8 +23,8 @@ function solve(model::CmeModel,
         matvec!(du, t, A, u)
         nothing
     end
-    fspprob = ODEProblem(odefun!, u0, tspan, p = model.parameters)
-    solutions = DE.solve(fspprob, ode_method, atol = odeatol, rtol = odertol, saveat = saveat)
+    fspprob = ODEProblem(odefun!, u0, tspan, model.parameters)
+    solutions = DE.solve(fspprob, ode_method, abstol = odeatol, reltol = odertol, saveat = saveat)
 
     output = FspOutputSparse{NS,IntT,RealT}(
         t = Vector{RealT}(),
@@ -106,12 +106,17 @@ function solve(model::CmeModel,
     initial_distribution::FspVectorSparse{NS,IntT,RealT},
     tspan::Tuple{AbstractFloat,AbstractFloat},
     fspalgorithm::AdaptiveFspSparse; saveat = [], fsptol::AbstractFloat = 1.0E-6,
-    odeatol::AbstractFloat = 1.0E-10,
+    odeatol::AbstractFloat = 1.0E-6,
     odertol::AbstractFloat = 1.0E-4,
     verbose::Bool = false) where {NS,IntT<:Integer,RealT<:AbstractFloat}
 
     tstart = min(tspan...)
     tend = max(tspan...)
+    if typeof(saveat) <: Number
+        _saveat = tspan[1]:saveat:tspan[2]
+    else
+        _saveat = saveat
+    end
     adapter = fspalgorithm.space_adapter
 
     p0 = deepcopy(initial_distribution.values)
@@ -150,8 +155,8 @@ function solve(model::CmeModel,
             abstol = eps()        
             )
 
-        fspprob = DE.ODEProblem(fsprhs!, unow, (tnow, tend), p = model.parameters)
-        integrator = DE.init(fspprob, fspalgorithm.ode_method, atol = odeatol, rtol = odertol, callback = fsp_cb, saveat = saveat)
+        fspprob = DE.ODEProblem(fsprhs!, unow, (tnow, tend), model.parameters)
+        integrator = DE.init(fspprob, fspalgorithm.ode_method, abstol = odeatol, reltol = odertol, callback = fsp_cb, saveat = _saveat)
 
         DE.step!(integrator, tend - tnow, true)
 

diff --git a/src/transientcme/sparse/rstepadapters.jl b/src/transientcme/sparse/rstepadapters.jl
@@ -98,7 +98,8 @@ function adapt!(statespace::StateSpaceSparse, adapter::SelectiveRStepAdapter,
     end
     sink_count = get_sink_count(statespace)
     du = similar(integrator.u)
-    integrator.f(du, integrator.u, [], t)
+    f! = integrator.f
+    f!(du, integrator.u, [], t)
     dsinks = du[end-sink_count+1:end]
     expandreactions = findall(dsinks .> 0)
     nold = get_state_count(statespace)