Merge pull request #438 from isaacsas/float-int-fix

fix #436

Author: ChrisRackauckas

src/systems/jumps/jumpsystem.jl

@@ -91,7 +91,7 @@ function assemble_crj(js, crj, statetoid)
 end
 
 function assemble_maj(js, maj::MassActionJump{U,Vector{Pair{V,W}},Vector{Pair{V2,W2}}},
-                      statetoid, subber) where {U,V,W,V2,W2}
+                      statetoid, subber, invttype) where {U,V,W,V2,W2}
     sr = maj.scaled_rates
     if sr isa Operation
         pval = subber(sr).value
@@ -118,7 +118,8 @@ function assemble_maj(js, maj::MassActionJump{U,Vector{Pair{V,W}},Vector{Pair{V2
     end
     sort!(ns)
 
-    MassActionJump(pval, rs, ns, scale_rates = false)
+    maj = MassActionJump(convert(invttype, pval), rs, ns, scale_rates = false)
+    return maj
 end
 
 """
@@ -163,12 +164,16 @@ sol = solve(jprob, SSAStepper())
 """
 function DiffEqJump.JumpProblem(js::JumpSystem, prob, aggregator; kwargs...)
 
-    statetoid = Dict(convert(Variable,state) => i for (i,state) in enumerate(states(js)))
-    parammap  = map((x,y)->Pair(x(),y), parameters(js), prob.p)
+    statetoid = Dict(convert(Variable,state) => i for (i,state) in enumerate(states(js)))    
     eqs       = equations(js)
-
-    subber = substituter(first.(parammap), last.(parammap))
-    majs = MassActionJump[assemble_maj(js, j, statetoid, subber) for j in eqs.x[1]]
+    invttype  = typeof(1 / prob.tspan[2])
+
+    # handling parameter substition and empty param vecs
+    p = (prob.p == DiffEqBase.NullParameters()) ? Operation[] : prob.p
+    parammap  = map((x,y)->Pair(x(),y), parameters(js), p)
+    subber    = substituter(first.(parammap), last.(parammap))
+    
+    majs = MassActionJump[assemble_maj(js, j, statetoid, subber, invttype) for j in eqs.x[1]]
     crjs = ConstantRateJump[assemble_crj(js, j, statetoid) for j in eqs.x[2]]
     vrjs = VariableRateJump[assemble_vrj(js, j, statetoid) for j in eqs.x[3]]
     ((prob isa DiscreteProblem) && !isempty(vrjs)) && error("Use continuous problems such as an ODEProblem or a SDEProblem with VariableRateJumps") 

src/systems/reaction/reactionsystem.jl

@@ -138,8 +138,10 @@ end
 function ReactionSystem(eqs, iv, species, params; systems = ReactionSystem[],
                                                   name = gensym(:ReactionSystem))
 
-    ReactionSystem(eqs, iv, convert.(Variable,species), convert.(Variable,params),
-                   name, systems)
+
+    isempty(species) && error("ReactionSystems require at least one species.")
+    paramvars = isempty(params) ? Variable[] : convert.(Variable, params)
+    ReactionSystem(eqs, iv, convert.(Variable,species), paramvars, name, systems)
 end
 
 # Calculate the ODE rate law
@@ -233,14 +235,19 @@ explicitly on the independent variable (usually time).
 """
 function ismassaction(rx, rs; rxvars = get_variables(rx.rate),
                               haveivdep = any(var -> isequal(rs.iv,convert(Variable,var)), rxvars))
-    return !(haveivdep || rx.only_use_rate || any(convert(Variable,rxv) in states(rs) for rxv in rxvars))
+    # if no dependencies must be zero order
+    if isempty(rxvars) 
+        return true
+    else
+        return !(haveivdep || rx.only_use_rate || any(convert(Variable,rxv) in states(rs) for rxv in rxvars))
+    end
 end
 
 function assemble_jumps(rs)
     eqs = Vector{Union{ConstantRateJump, MassActionJump, VariableRateJump}}()
 
-    for rx in equations(rs)
-        rxvars    = get_variables(rx.rate)
+    for rx in equations(rs)        
+        rxvars    = (rx.rate isa Operation) ? get_variables(rx.rate) : Operation[]
         haveivdep = any(var -> isequal(rs.iv,convert(Variable,var)), rxvars)
         if ismassaction(rx, rs; rxvars=rxvars, haveivdep=haveivdep)
             reactant_stoch = isempty(rx.substoich) ? [0 => 1] : [var2op(sub.op) => stoich for (sub,stoich) in zip(rx.substrates,rx.substoich)]

test/reactionsystem.jl

@@ -88,7 +88,7 @@ js = convert(JumpSystem, rs)
 @test all(map(i -> typeof(js.eqs[i]) <: DiffEqJump.ConstantRateJump, 15:18))
 @test all(map(i -> typeof(js.eqs[i]) <: DiffEqJump.VariableRateJump, 19:20))
 
-pars = rand(length(k)); u0 = rand(1:100,4); time = rand();
+pars = rand(length(k)); u0 = rand(1:10,4); time = rand();
 jumps = Vector{Union{ConstantRateJump, MassActionJump, VariableRateJump}}(undef,length(js.eqs))
 
 jumps[1] = MassActionJump(pars[1], Vector{Pair{Int,Int}}(), [1 => 1]);
@@ -117,23 +117,36 @@ jumps[20] = VariableRateJump((u,p,t) -> p[20]*t*u[1]*binomial(u[2],2)*u[3], inte
 statetoid = Dict(convert(Variable,state) => i for (i,state) in enumerate(states(js)))
 parammap = map((x,y)->Pair(x(),y),parameters(js),pars)
 for i = 1:14
-  maj = MT.assemble_maj(js, js.eqs[i], statetoid, ModelingToolkit.substituter(first.(parammap), last.(parammap)))
+  maj = MT.assemble_maj(js, js.eqs[i], statetoid, ModelingToolkit.substituter(first.(parammap), last.(parammap)),eltype(pars))
   @test abs(jumps[i].scaled_rates - maj.scaled_rates) < 100*eps()
   @test jumps[i].reactant_stoch == maj.reactant_stoch
   @test jumps[i].net_stoch == maj.net_stoch
 end
 for i = 15:18
   (i==16) && continue
   crj = MT.assemble_crj(js, js.eqs[i], statetoid)
-  @test abs(crj.rate(u0,p,time) - jumps[i].rate(u0,p,time)) < 100*eps()
+  @test isapprox(crj.rate(u0,p,time), jumps[i].rate(u0,p,time))
   fake_integrator1 = (u=zeros(4),p=p,t=0); fake_integrator2 = deepcopy(fake_integrator1);
   crj.affect!(fake_integrator1); jumps[i].affect!(fake_integrator2);
   @test fake_integrator1 == fake_integrator2
 end
 for i = 19:20
   crj = MT.assemble_vrj(js, js.eqs[i], statetoid)
-  @test abs(crj.rate(u0,p,time) - jumps[i].rate(u0,p,time)) < 100*eps()
+  @test isapprox(crj.rate(u0,p,time), jumps[i].rate(u0,p,time))
   fake_integrator1 = (u=zeros(4),p=p,t=0.); fake_integrator2 = deepcopy(fake_integrator1);
   crj.affect!(fake_integrator1); jumps[i].affect!(fake_integrator2);
   @test fake_integrator1 == fake_integrator2
 end
+
+
+# test for https://github.com/SciML/ModelingToolkit.jl/issues/436
+@parameters t 
+@variables S I
+rxs = [Reaction(1,[S],[I]), Reaction(1.1,[S],[I])]
+rs = ReactionSystem(rxs, t, [S,I], [])
+js = convert(JumpSystem, rs)
+dprob = DiscreteProblem(js, [S => 1, I => 1], (0.0,10.0))
+jprob = JumpProblem(js, dprob, Direct())
+sol = solve(jprob, SSAStepper())
+
+nothing