From 7a95021fc87f8c303d79e4a166d40d0c720d5e82 Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 16 Oct 2024 16:56:05 -0400 Subject: [PATCH] format --- docs/make.jl | 34 ++--- docs/old/getting_started.md | 2 +- docs/pages.jl | 12 +- .../applications/advanced_point_process.md | 112 ++++++++-------- docs/src/faq.md | 6 +- docs/src/jump_types.md | 10 +- .../tutorials/discrete_stochastic_example.md | 4 +- .../src/tutorials/point_process_simulation.md | 2 +- docs/src/tutorials/spatial.md | 6 +- src/SSA_stepper.jl | 46 +++---- src/aggregators/aggregated_api.jl | 16 +-- src/aggregators/aggregators.jl | 8 +- src/aggregators/bracketing.jl | 2 +- src/aggregators/ccnrm.jl | 24 ++-- src/aggregators/coevolve.jl | 34 ++--- src/aggregators/direct.jl | 20 +-- src/aggregators/directcr.jl | 26 ++-- src/aggregators/frm.jl | 24 ++-- src/aggregators/nrm.jl | 20 +-- src/aggregators/prioritytable.jl | 6 +- src/aggregators/rdirect.jl | 27 ++-- src/aggregators/rssa.jl | 24 ++-- src/aggregators/rssacr.jl | 30 ++--- src/aggregators/sortingdirect.jl | 20 +-- src/aggregators/ssajump.jl | 27 ++-- src/coupling.jl | 20 +-- src/extended_jump_array.jl | 54 ++++---- src/jumps.jl | 126 +++++++++--------- src/massaction_rates.jl | 20 +-- src/problem.jl | 104 +++++++-------- src/simple_regular_solve.jl | 8 +- src/solve.jl | 28 ++-- src/spatial/bracketing.jl | 2 +- src/spatial/directcrdirect.jl | 48 +++---- src/spatial/flatten.jl | 36 ++--- src/spatial/hop_rates.jl | 85 ++++++------ src/spatial/nsm.jl | 47 +++---- src/spatial/reaction_rates.jl | 6 +- src/spatial/spatial_massaction_jump.jl | 80 +++++------ src/spatial/topology.jl | 8 +- src/spatial/utils.jl | 6 +- test/allocations.jl | 12 +- test/bimolerx_test.jl | 20 +-- test/bracketing.jl | 2 +- test/degenerate_rx_cases.jl | 10 +- test/extended_jump_array.jl | 6 +- test/extinction_test.jl | 6 +- test/functionwrappers.jl | 2 +- test/geneexpr_test.jl | 40 +++--- test/hawkes_test.jl | 12 +- test/jprob_symbol_indexing.jl | 2 +- test/linearreaction_test.jl | 22 +-- test/monte_carlo_test.jl | 4 +- test/remake_test.jl | 2 +- test/reversible_binding.jl | 6 +- test/save_positions.jl | 8 +- test/saveat_regression.jl | 2 +- test/spatial/ABC.jl | 14 +- test/spatial/bracketing.jl | 2 +- test/spatial/diffusion.jl | 70 +++++----- test/spatial/hop_rates.jl | 26 ++-- test/spatial/spatial_majump.jl | 42 +++--- test/spatial/topology.jl | 2 +- test/splitcoupled.jl | 27 ++-- test/ssa_callback_test.jl | 22 +-- test/ssa_tests.jl | 2 +- test/thread_safety.jl | 2 +- test/variable_rate.jl | 16 +-- 68 files changed, 807 insertions(+), 794 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index c362e135..46825545 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -8,24 +8,24 @@ cp(joinpath(docpath, "Project.toml"), joinpath(assetpath, "Project.toml"), force include("pages.jl") mathengine = MathJax3(Dict(:loader => Dict("load" => ["[tex]/require", "[tex]/mathtools"]), - :tex => Dict("inlineMath" => [["\$", "\$"], ["\\(", "\\)"]], - "packages" => [ - "base", - "ams", - "autoload", - "mathtools", - "require", - ]))) + :tex => Dict("inlineMath" => [["\$", "\$"], ["\\(", "\\)"]], + "packages" => [ + "base", + "ams", + "autoload", + "mathtools", + "require" + ]))) makedocs(sitename = "JumpProcesses.jl", - authors = "Chris Rackauckas", - modules = [JumpProcesses], - clean = true, doctest = false, linkcheck = true, warnonly = [:missing_docs], - format = Documenter.HTML(; assets = ["assets/favicon.ico"], - canonical = "https://docs.sciml.ai/JumpProcesses/", - prettyurls = (get(ENV, "CI", nothing) == "true"), - mathengine), - pages = pages) + authors = "Chris Rackauckas", + modules = [JumpProcesses], + clean = true, doctest = false, linkcheck = true, warnonly = [:missing_docs], + format = Documenter.HTML(; assets = ["assets/favicon.ico"], + canonical = "https://docs.sciml.ai/JumpProcesses/", + prettyurls = (get(ENV, "CI", nothing) == "true"), + mathengine), + pages = pages) deploydocs(repo = "github.com/SciML/JumpProcesses.jl.git"; - push_preview = true) + push_preview = true) diff --git a/docs/old/getting_started.md b/docs/old/getting_started.md index 7f6534bf..3192de0c 100644 --- a/docs/old/getting_started.md +++ b/docs/old/getting_started.md @@ -44,7 +44,7 @@ sol = solve(jprob, SSAStepper()) using Plots plot(sol, title = "Sample path from a jump process with constant rate", - label = "N(t)", xlabel = "t", legend = :bottomright) + label = "N(t)", xlabel = "t", legend = :bottomright) ``` ## Step 1: Defining a problem diff --git a/docs/pages.jl b/docs/pages.jl index 5d79940d..1f973f24 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -2,13 +2,13 @@ pages = ["index.md", "Tutorials" => Any["tutorials/simple_poisson_process.md", - "tutorials/discrete_stochastic_example.md", - "tutorials/point_process_simulation.md", - "tutorials/jump_diffusion.md", - "tutorials/spatial.md"], + "tutorials/discrete_stochastic_example.md", + "tutorials/point_process_simulation.md", + "tutorials/jump_diffusion.md", + "tutorials/spatial.md"], "Applications" => Any["applications/advanced_point_process.md"], "Type Documentation" => Any["Jumps, JumpProblem, and Aggregators" => "jump_types.md", - "Jump solvers" => "jump_solve.md"], + "Jump solvers" => "jump_solve.md"], "FAQ" => "faq.md", - "API" => "api.md", + "API" => "api.md" ] diff --git a/docs/src/applications/advanced_point_process.md b/docs/src/applications/advanced_point_process.md index 3e122923..1b7ab93c 100644 --- a/docs/src/applications/advanced_point_process.md +++ b/docs/src/applications/advanced_point_process.md @@ -91,7 +91,7 @@ end function Base.show(io::IO, pp::SciMLPointProcess) println(io, - "SciMLPointProcess with $(length(pp.jumps)) processes on the interval [$(pp.tmin), $(pp.tmax)).") + "SciMLPointProcess with $(length(pp.jumps)) processes on the interval [$(pp.tmin), $(pp.tmax)).") end ``` @@ -294,12 +294,12 @@ function hawkes_p(pp::SciMLPointProcess{M, J, G, D, T}) where {M, J, G, D, T} tmax = pp.tmax h = History(; times = T[], marks = Tuple{Int, M}[], tmin = tmin, tmax = tmax) return (λ = 0.5, - α = 0.1, - β = 2.0, - ϕ = zeros(T, length(g)), - M = map(d -> tuple(mean(d)...), pp.mark_dist), - T = zeros(T, length(g)), - h = h) + α = 0.1, + β = 2.0, + ϕ = zeros(T, length(g)), + M = map(d -> tuple(mean(d)...), pp.mark_dist), + T = zeros(T, length(g)), + h = h) end function hawkes_p(pp::SciMLPointProcess{M, J, G, D, T}, p) where {M, J, G, D, T} @@ -324,17 +324,17 @@ mark_dist = [MvNormal(rand(2), [0.2, 0.2]) for i in 1:nv(G)] jumps = [hawkes_jump(i, g, mark_dist) for i in 1:nv(G)] tspan = (0.0, 50.0) hawkes = SciMLPointProcess{ - Vector{Real}, - eltype(jumps), - typeof(g), - eltype(mark_dist), - eltype(tspan) - }(jumps, - mark_dist, - g, - hawkes_p, - tspan[1], - tspan[2]) + Vector{Real}, + eltype(jumps), + typeof(g), + eltype(mark_dist), + eltype(tspan) +}(jumps, + mark_dist, + g, + hawkes_p, + tspan[1], + tspan[2]) ``` ## [Sampling](@id tpp_sampling) @@ -378,10 +378,10 @@ function Base.rand(pp::SciMLPointProcess, n::Int) end function Base.rand(rng::AbstractRNG, - pp::SciMLPointProcess{M, J, G, D, T}, - tmin::T, - tmax::T, - n::Int) where {M, J, G, D, T <: Real} + pp::SciMLPointProcess{M, J, G, D, T}, + tmin::T, + tmax::T, + n::Int) where {M, J, G, D, T <: Real} tspan = (tmin, tmax) save_positions = (false, false) out = Array{History, 1}(undef, n) @@ -487,7 +487,7 @@ function intensity(pp::SciMLPointProcess, t, h; saveat = [], save_positions = (t end callback = DiscreteCallback(condition, affect!; save_positions) dprob = DiscreteProblem(rates, rates(nothing, p, tmin), (tmin, tmax), p; callback, - tstops, saveat) + tstops, saveat) sol = solve(dprob, FunctionMap()) return sol end @@ -573,11 +573,11 @@ conditional intensity using an ODEProblem. ```@example tpp-advanced function integrated_intensity(pp::SciMLPointProcess, - t, - h; - alg = nothing, - saveat = [], - save_positions = (true, true)) + t, + h; + alg = nothing, + saveat = [], + save_positions = (true, true)) p = params(pp) times = event_times(h) marks = event_marks(h) @@ -593,12 +593,12 @@ function integrated_intensity(pp::SciMLPointProcess, end callback = DiscreteCallback(condition, affect!; save_positions) prob = ODEProblem(rates, - zeros(eltype(times), length(pp.jumps)), - tspan, - p; - tstops = times, - callback, - saveat) + zeros(eltype(times), length(pp.jumps)), + tspan, + p; + tstops = times, + callback, + saveat) sol = solve(prob, alg) return sol end @@ -613,11 +613,11 @@ the next section when we compute the log-likelihood. ```@example tpp-advanced function integrated_ground_intensity(pp::SciMLPointProcess, h, a, b) Λ = integrated_intensity(pp, - b, - h; - alg = Rodas4P(), - saveat = [a, b], - save_positions = (false, false)) + b, + h; + alg = Rodas4P(), + saveat = [a, b], + save_positions = (false, false)) return sum(Λ(b)) - sum(Λ(a)) end @@ -631,10 +631,10 @@ solver that can deal with stiff problems like `Rodas4P()`. ```@example tpp-advanced Λ = integrated_intensity(hawkes, - max_time(hawkes), - h; - saveat = event_times(h), - alg = Rodas4P()) + max_time(hawkes), + h; + saveat = event_times(h), + alg = Rodas4P()) plot(Λ) ``` @@ -662,7 +662,7 @@ function hawkes_integrated_intensity(pp::SciMLPointProcess, t, h; saveat = []) return u end dprob = DiscreteProblem(compensator, zeros(eltype(event_times(h)), length(pp.jumps)), - (min_time(h), t), p; tstops, saveat) + (min_time(h), t), p; tstops, saveat) sol = solve(dprob, FunctionMap()) return sol end @@ -693,10 +693,10 @@ function Base.filter(f, h::History) end end return History(; - times = filtered_times, - marks = filtered_marks, - tmin = min_time(h), - tmax = max_time(h)) + times = filtered_times, + marks = filtered_marks, + tmin = min_time(h), + tmax = max_time(h)) end ``` @@ -744,10 +744,10 @@ itself to produce the QQ-plot for the ground process. for _ in 1:250 _h = rand(hawkes) _Λ = integrated_intensity(hawkes, - max_time(hawkes), - _h; - saveat = event_times(_h), - alg = Rodas4P()) + max_time(hawkes), + _h; + saveat = event_times(_h), + alg = Rodas4P()) _h̃ = time_change(_h, (t) -> sum(_Λ(t))) append!(Δt̃, diff(event_times(_h̃))) end @@ -763,10 +763,10 @@ Likewise, we can produce the QQ-plot for each sub-TPP. for _ in 1:250 _h = rand(hawkes) _Λ = integrated_intensity(hawkes, - max_time(hawkes), - _h; - saveat = event_times(_h), - alg = Rodas4P()) + max_time(hawkes), + _h; + saveat = event_times(_h), + alg = Rodas4P()) for i in 1:nv(G) _h̃ = time_change(filter((t, mark) -> mark[1] == i, _h), (t) -> Λ(t; idxs = i)) append!(Δt̃[i], diff(event_times(_h̃))) diff --git a/docs/src/faq.md b/docs/src/faq.md index 1875d773..d2d0be3b 100644 --- a/docs/src/faq.md +++ b/docs/src/faq.md @@ -50,7 +50,7 @@ vj2 = VariableRateJump(rate4, affect4!) vjtuple = (vj1, vj2) jset = JumpSet(; constant_jumps = cjvec, variable_jumps = vjtuple, - massaction_jumps = mass_act_jump) + massaction_jumps = mass_act_jump) ``` ## How can I set the random number generator used in the jump process sampling algorithms (SSAs)? @@ -62,7 +62,7 @@ argument. Continuing the previous example: #] add RandomNumbers using RandomNumbers jprob = JumpProblem(dprob, Direct(), maj, - rng = Xorshifts.Xoroshiro128Star(rand(UInt64))) + rng = Xorshifts.Xoroshiro128Star(rand(UInt64))) ``` uses the `Xoroshiro128Star` generator from @@ -148,7 +148,7 @@ function paffect!(integrator) nothing end sol = solve(jprob, SSAStepper(), tstops = [20.0], - callback = DiscreteCallback(pcondit, paffect!)) + callback = DiscreteCallback(pcondit, paffect!)) ``` Here at time `20.0` we turn off production of `u[2]` while activating production diff --git a/docs/src/jump_types.md b/docs/src/jump_types.md index faae6230..2d0e6481 100644 --- a/docs/src/jump_types.md +++ b/docs/src/jump_types.md @@ -161,9 +161,9 @@ The constructor for a [`VariableRateJump`](@ref) is: ```julia VariableRateJump(rate, affect!; - lrate = nothing, urate = nothing, rateinterval = nothing, - idxs = nothing, rootfind = true, save_positions = (true, true), - interp_points = 10, abstol = 1e-12, reltol = 0) + lrate = nothing, urate = nothing, rateinterval = nothing, + idxs = nothing, rootfind = true, save_positions = (true, true), + interp_points = 10, abstol = 1e-12, reltol = 0) ``` - `rate(u, p, t)` is a function which calculates the rate given the current @@ -241,8 +241,8 @@ is: ```julia JumpProblem(prob, aggregator, jumps::JumpSet; - save_positions = prob isa AbstractDiscreteProblem ? (false, true) : - (true, true)) + save_positions = prob isa AbstractDiscreteProblem ? (false, true) : + (true, true)) ``` The aggregator is the method for simulating `ConstantRateJump`s, diff --git a/docs/src/tutorials/discrete_stochastic_example.md b/docs/src/tutorials/discrete_stochastic_example.md index 73724aaf..3638e70e 100644 --- a/docs/src/tutorials/discrete_stochastic_example.md +++ b/docs/src/tutorials/discrete_stochastic_example.md @@ -453,11 +453,11 @@ latter determining the affect function. rateidxs = [1, 2] # i.e., [β, ν] reactant_stoich = [ [1 => 1, 2 => 1], # 1*S and 1*I - [2 => 1], # 1*I + [2 => 1] # 1*I ] net_stoich = [ [1 => -1, 2 => 1], # -1*S and 1*I - [2 => -1, 3 => 1], # -1*I and 1*R + [2 => -1, 3 => 1] # -1*I and 1*R ] mass_act_jump = MassActionJump(reactant_stoich, net_stoich; param_idxs = rateidxs) ``` diff --git a/docs/src/tutorials/point_process_simulation.md b/docs/src/tutorials/point_process_simulation.md index 83d85e25..40e6d5f5 100644 --- a/docs/src/tutorials/point_process_simulation.md +++ b/docs/src/tutorials/point_process_simulation.md @@ -196,7 +196,7 @@ Now, we can declare the seasonal process as a `VariableRateJump`. ```@example tpp-tutorial seasonal_process = VariableRateJump(seasonal_rate, seasonal_affect!; - urate, rateinterval, lrate) + urate, rateinterval, lrate) ``` We initialize a new base problem with a different simulation diff --git a/docs/src/tutorials/spatial.md b/docs/src/tutorials/spatial.md index 2b48257a..314b2e47 100644 --- a/docs/src/tutorials/spatial.md +++ b/docs/src/tutorials/spatial.md @@ -73,7 +73,7 @@ We are now ready to set up the `JumpProblem` with the Next Subvolume Method. ```@example spatial alg = NSM() jump_prob = JumpProblem(prob, alg, majumps, hopping_constants = hopping_constants, - spatial_system = grid, save_positions = (true, false)) + spatial_system = grid, save_positions = (true, false)) ``` The `save_positions` keyword tells the solver to save the positions just before the jumps. To solve the jump problem do @@ -124,7 +124,7 @@ function get_frame(k, sol, linear_size, labels, title) end for species in 1:num_species scatter!(plt, species_seriess_x[species], species_seriess_y[species], - label = labels[species], marker = 6) + label = labels[species], marker = 6) end xticks!(plt, range(xlim..., length = linear_size + 1)) yticks!(plt, range(ylim..., length = linear_size + 1)) @@ -146,7 +146,7 @@ function animate_2d(sol, linear_size; species_labels, title, verbose = true) end # animate anim = animate_2d(solution, 5, species_labels = ["A", "B", "C"], title = "A + B <--> C", - verbose = false) + verbose = false) fps = 5 gif(anim, fps = fps) ``` diff --git a/src/SSA_stepper.jl b/src/SSA_stepper.jl index da0c9e84..82c3c93d 100644 --- a/src/SSA_stepper.jl +++ b/src/SSA_stepper.jl @@ -135,7 +135,7 @@ function DiffEqBase.solve!(integrator::SSAIntegrator) # check callbacks one last time if !(integrator.opts.callback.discrete_callbacks isa Tuple{}) DiffEqBase.apply_discrete_callback!(integrator, - integrator.opts.callback.discrete_callbacks...) + integrator.opts.callback.discrete_callbacks...) end if integrator.saveat !== nothing && !isempty(integrator.saveat) @@ -162,15 +162,15 @@ function DiffEqBase.solve!(integrator::SSAIntegrator) end function DiffEqBase.__init(jump_prob::JumpProblem, - alg::SSAStepper; - save_start = true, - save_end = true, - seed = nothing, - alias_jump = Threads.threadid() == 1, - saveat = nothing, - callback = nothing, - tstops = nothing, - numsteps_hint = 100) + alg::SSAStepper; + save_start = true, + save_end = true, + seed = nothing, + alias_jump = Threads.threadid() == 1, + saveat = nothing, + callback = nothing, + tstops = nothing, + numsteps_hint = 100) # hack until alias system is in place alias_tstops = false @@ -207,9 +207,9 @@ function DiffEqBase.__init(jump_prob::JumpProblem, save_everystep = any(cb.save_positions) sol = DiffEqBase.build_solution(prob, alg, t, u, dense = save_everystep, - calculate_error = false, - stats = DiffEqBase.Stats(0), - interp = DiffEqBase.ConstantInterpolation(t, u)) + calculate_error = false, + stats = DiffEqBase.Stats(0), + interp = DiffEqBase.ConstantInterpolation(t, u)) _saveat = (saveat isa Number) ? (prob.tspan[1]:saveat:prob.tspan[2]) : saveat if _saveat !== nothing && !isempty(_saveat) && _saveat[1] == prob.tspan[1] @@ -241,9 +241,9 @@ function DiffEqBase.__init(jump_prob::JumpProblem, end integrator = SSAIntegrator(prob.f, copy(prob.u0), prob.tspan[1], prob.tspan[1], tdir, - prob.p, sol, 1, prob.tspan[1], cb, _saveat, save_everystep, - save_end, cur_saveat, opts, _tstops, 1, false, true, - alias_tstops, false) + prob.p, sol, 1, prob.tspan[1], cb, _saveat, save_everystep, + save_end, cur_saveat, opts, _tstops, 1, false, true, + alias_tstops, false) cb.initialize(cb, integrator.u, prob.tspan[1], integrator) DiffEqBase.initialize!(opts.callback, integrator.u, prob.tspan[1], integrator) integrator @@ -285,7 +285,7 @@ end # The Jump aggregators should not register the next jump through add_tstop! for SSAIntegrator # such that we can achieve maximum performance @inline function register_next_jump_time!(integrator::SSAIntegrator, - p::AbstractSSAJumpAggregator, t) + p::AbstractSSAJumpAggregator, t) integrator.tstop = p.next_jump_time nothing end @@ -330,7 +330,7 @@ function DiffEqBase.step!(integrator::SSAIntegrator) if !(integrator.opts.callback.discrete_callbacks isa Tuple{}) discrete_modified, saved_in_cb = DiffEqBase.apply_discrete_callback!(integrator, - integrator.opts.callback.discrete_callbacks...) + integrator.opts.callback.discrete_callbacks...) else saved_in_cb = false end @@ -384,10 +384,10 @@ end export SSAStepper function SciMLBase.isdenseplot(sol::ODESolution{ - T, N, uType, uType2, DType, tType, rateType, - discType, P, - SSAStepper}) where {T, N, uType, uType2, - DType, tType, rateType, - discType, P} + T, N, uType, uType2, DType, tType, rateType, + discType, P, + SSAStepper}) where {T, N, uType, uType2, + DType, tType, rateType, + discType, P} sol.dense end diff --git a/src/aggregators/aggregated_api.jl b/src/aggregators/aggregated_api.jl index 82f1f7fc..a65c4ad8 100644 --- a/src/aggregators/aggregated_api.jl +++ b/src/aggregators/aggregated_api.jl @@ -11,28 +11,28 @@ Notes vector is unchanged, this can safely be set to false to improve performance. """ function reset_aggregated_jumps!(integrator, uprev = nothing; update_jump_params = true, - kwargs...) + kwargs...) reset_aggregated_jumps!(integrator, uprev, integrator.opts.callback, - update_jump_params = update_jump_params, kwargs...) + update_jump_params = update_jump_params, kwargs...) nothing end function reset_aggregated_jumps!(integrator, uprev, callback::Nothing; - update_jump_params = true, kwargs...) + update_jump_params = true, kwargs...) nothing end function reset_aggregated_jumps!(integrator, uprev, callback::CallbackSet; - update_jump_params = true, kwargs...) + update_jump_params = true, kwargs...) if !isempty(callback.discrete_callbacks) reset_aggregated_jumps!(integrator, uprev, callback.discrete_callbacks..., - update_jump_params = update_jump_params, kwargs...) + update_jump_params = update_jump_params, kwargs...) end nothing end function reset_aggregated_jumps!(integrator, uprev, cb::DiscreteCallback, cbs...; - update_jump_params = true, kwargs...) + update_jump_params = true, kwargs...) if cb.condition isa AbstractSSAJumpAggregator maj = cb.condition.ma_jumps update_jump_params && using_params(maj) && @@ -40,12 +40,12 @@ function reset_aggregated_jumps!(integrator, uprev, cb::DiscreteCallback, cbs... cb.condition(cb, integrator.u, integrator.t, integrator) end reset_aggregated_jumps!(integrator, uprev, cbs...; - update_jump_params = update_jump_params, kwargs...) + update_jump_params = update_jump_params, kwargs...) nothing end function reset_aggregated_jumps!(integrator, uprev, cb::DiscreteCallback; - update_jump_params = true, kwargs...) + update_jump_params = true, kwargs...) if cb.condition isa AbstractSSAJumpAggregator maj = cb.condition.ma_jumps update_jump_params && using_params(maj) && diff --git a/src/aggregators/aggregators.jl b/src/aggregators/aggregators.jl index 4362c07e..732a542c 100644 --- a/src/aggregators/aggregators.jl +++ b/src/aggregators/aggregators.jl @@ -173,7 +173,7 @@ algorithm with optimal binning, Journal of Chemical Physics 143, 074108 struct DirectCRDirect <: AbstractAggregatorAlgorithm end const JUMP_AGGREGATORS = (Direct(), DirectFW(), DirectCR(), SortingDirect(), RSSA(), FRM(), - FRMFW(), NRM(), RSSACR(), RDirect(), Coevolve(), CCNRM()) + FRMFW(), NRM(), RSSACR(), RDirect(), Coevolve(), CCNRM()) # For JumpProblem construction without an aggregator struct NullAggregator <: AbstractAggregatorAlgorithm end @@ -205,9 +205,9 @@ is_spatial(aggregator::DirectCRDirect) = true # return the fastest aggregator out of the available ones function select_aggregator(jumps::JumpSet; vartojumps_map = nothing, - jumptovars_map = nothing, dep_graph = nothing, - spatial_system = nothing, - hopping_constants = nothing) + jumptovars_map = nothing, dep_graph = nothing, + spatial_system = nothing, + hopping_constants = nothing) # detect if a spatial SSA should be used !isnothing(spatial_system) && !isnothing(hopping_constants) && return DirectCRDirect diff --git a/src/aggregators/bracketing.jl b/src/aggregators/bracketing.jl index ba0c0712..40146c7a 100644 --- a/src/aggregators/bracketing.jl +++ b/src/aggregators/bracketing.jl @@ -69,7 +69,7 @@ get brackets for the rate of reaction rx by first checking if the reaction is a return get_majump_brackets(p.ulow, p.uhigh, rx, ma_jumps) else @inbounds return get_cjump_brackets(p.ulow, p.uhigh, p.rates[rx - num_majumps], - params, t) + params, t) end end diff --git a/src/aggregators/ccnrm.jl b/src/aggregators/ccnrm.jl index 5a244512..3adfb2ed 100644 --- a/src/aggregators/ccnrm.jl +++ b/src/aggregators/ccnrm.jl @@ -21,9 +21,9 @@ mutable struct CCNRMJumpAggregation{T, S, F1, F2, RNG, DEPGR, PT} <: end function CCNRMJumpAggregation(nj::Int, njt::T, et::T, crs::Vector{T}, sr::T, - maj::S, rs::F1, affs!::F2, sps::Tuple{Bool, Bool}, - rng::RNG; num_specs, dep_graph = nothing, - kwargs...) where {T, S, F1, F2, RNG} + maj::S, rs::F1, affs!::F2, sps::Tuple{Bool, Bool}, + rng::RNG; num_specs, dep_graph = nothing, + kwargs...) where {T, S, F1, F2, RNG} # a dependency graph is needed and must be provided if there are constant rate jumps if dep_graph === nothing @@ -42,27 +42,27 @@ function CCNRMJumpAggregation(nj::Int, njt::T, et::T, crs::Vector{T}, sr::T, binwidthconst = haskey(kwargs, :binwidthconst) ? kwargs[:binwidthconst] : 16 numbinsconst = haskey(kwargs, :numbinsconst) ? kwargs[:numbinsconst] : 20 ptt = PriorityTimeTable(zeros(T, length(crs)), zero(T), one(T), - binwidthconst = binwidthconst, numbinsconst = numbinsconst) # We will re-initialize this in initialize!() + binwidthconst = binwidthconst, numbinsconst = numbinsconst) # We will re-initialize this in initialize!() affecttype = F2 <: Tuple ? F2 : Any CCNRMJumpAggregation{T, S, F1, affecttype, RNG, typeof(dg), typeof(ptt)}(nj, nj, njt, - et, - crs, sr, maj, - rs, affs!, sps, - rng, dg, ptt) + et, + crs, sr, maj, + rs, affs!, sps, + rng, dg, ptt) end +############################# Required Functions ############################## # creating the JumpAggregation structure (function wrapper-based constant jumps) function aggregate(aggregator::CCNRM, u, p, t, end_time, constant_jumps, - ma_jumps, save_positions, rng; kwargs...) + ma_jumps, save_positions, rng; kwargs...) # handle constant jumps using function wrappers rates, affects! = get_jump_info_fwrappers(u, p, t, constant_jumps) build_jump_aggregation(CCNRMJumpAggregation, u, p, t, end_time, ma_jumps, - rates, affects!, save_positions, rng; num_specs = length(u), - kwargs...) + rates, affects!, save_positions, rng; num_specs = length(u), + kwargs...) end # set up a new simulation and calculate the first jump / jump time @@ -114,7 +114,7 @@ function update_dependent_rates!(p::CCNRMJumpAggregation, u, params, t) # update the jump rate @inbounds cur_rates[rx] = calculate_jump_rate(ma_jumps, num_majumps, rates, u, - params, t, rx) + params, t, rx) # Calculate new jump times for dependent jumps if rx != p.next_jump && oldrate > zero(oldrate) diff --git a/src/aggregators/coevolve.jl b/src/aggregators/coevolve.jl index 78a8230f..82e8beb9 100644 --- a/src/aggregators/coevolve.jl +++ b/src/aggregators/coevolve.jl @@ -24,10 +24,10 @@ mutable struct CoevolveJumpAggregation{T, S, F1, F2, RNG, GR, PQ} <: end function CoevolveJumpAggregation(nj::Int, njt::T, et::T, crs::Vector{T}, sr::Nothing, - maj::S, rs::F1, affs!::F2, sps::Tuple{Bool, Bool}, - rng::RNG; u::U, dep_graph = nothing, lrates, urates, - rateintervals, haslratevec, - cur_lrates::Vector{T}) where {T, S, F1, F2, RNG, U} + maj::S, rs::F1, affs!::F2, sps::Tuple{Bool, Bool}, + rng::RNG; u::U, dep_graph = nothing, lrates, urates, + rateintervals, haslratevec, + cur_lrates::Vector{T}) where {T, S, F1, F2, RNG, U} if dep_graph === nothing if (get_num_majumps(maj) == 0) || !isempty(urates) error("To use Coevolve a dependency graph between jumps must be supplied.") @@ -50,10 +50,10 @@ function CoevolveJumpAggregation(nj::Int, njt::T, et::T, crs::Vector{T}, sr::Not pq = MutableBinaryMinHeap{T}() affecttype = F2 <: Tuple ? F2 : Any CoevolveJumpAggregation{T, S, F1, affecttype, RNG, typeof(dg), - typeof(pq)}(nj, nj, njt, et, crs, sr, maj, - rs, affs!, sps, rng, dg, pq, - lrates, urates, rateintervals, - haslratevec, cur_lrates) + typeof(pq)}(nj, nj, njt, et, crs, sr, maj, + rs, affs!, sps, rng, dg, pq, + lrates, urates, rateintervals, + haslratevec, cur_lrates) end # display @@ -85,8 +85,8 @@ function (p::CoevolveJumpAggregation)(integrator::I) where {I <: end function (p::CoevolveJumpAggregation{ - T, S, F1, F2})(integrator::AbstractSSAIntegrator) where - {T, S, F1, F2 <: Union{Tuple, Nothing}} + T, S, F1, F2})(integrator::AbstractSSAIntegrator) where + {T, S, F1, F2 <: Union{Tuple, Nothing}} if !accept_next_jump!(p, integrator, integrator.u, integrator.p, integrator.t) return nothing end @@ -98,10 +98,10 @@ end # creating the JumpAggregation structure (tuple-based variable jumps) function aggregate(aggregator::Coevolve, u, p, t, end_time, constant_jumps, - ma_jumps, save_positions, rng; dep_graph = nothing, - variable_jumps = nothing, kwargs...) + ma_jumps, save_positions, rng; dep_graph = nothing, + variable_jumps = nothing, kwargs...) RateWrapper = FunctionWrappers.FunctionWrapper{typeof(t), - Tuple{typeof(u), typeof(p), typeof(t)}} + Tuple{typeof(u), typeof(p), typeof(t)}} ncrjs = (constant_jumps === nothing) ? 0 : length(constant_jumps) nvrjs = (variable_jumps === nothing) ? 0 : length(variable_jumps) @@ -141,9 +141,9 @@ function aggregate(aggregator::Coevolve, u, p, t, end_time, constant_jumps, next_jump = 0 next_jump_time = typemax(t) CoevolveJumpAggregation(next_jump, next_jump_time, end_time, cur_rates, sum_rate, - ma_jumps, rates, affects!, save_positions, rng; - u, dep_graph, lrates, urates, rateintervals, haslratevec, - cur_lrates) + ma_jumps, rates, affects!, save_positions, rng; + u, dep_graph, lrates, urates, rateintervals, haslratevec, + cur_lrates) end # set up a new simulation and calculate the first jump / jump time @@ -156,7 +156,7 @@ end # execute one jump, changing the system state function execute_jumps!(p::CoevolveJumpAggregation, integrator, u, params, t, - affects!) + affects!) # execute jump update_state!(p, integrator, u, affects!) # update current jump rates and times diff --git a/src/aggregators/direct.jl b/src/aggregators/direct.jl index 028e5603..ab33dd84 100644 --- a/src/aggregators/direct.jl +++ b/src/aggregators/direct.jl @@ -13,35 +13,35 @@ mutable struct DirectJumpAggregation{T, S, F1, F2, RNG} <: rng::RNG end function DirectJumpAggregation(nj::Int, njt::T, et::T, crs::Vector{T}, sr::T, maj::S, - rs::F1, affs!::F2, sps::Tuple{Bool, Bool}, rng::RNG; - kwargs...) where {T, S, F1, F2, RNG} + rs::F1, affs!::F2, sps::Tuple{Bool, Bool}, rng::RNG; + kwargs...) where {T, S, F1, F2, RNG} affecttype = F2 <: Tuple ? F2 : Any DirectJumpAggregation{T, S, F1, affecttype, RNG}(nj, nj, njt, et, crs, sr, maj, rs, - affs!, sps, rng) + affs!, sps, rng) end ############################# Required Functions ############################# # creating the JumpAggregation structure (tuple-based constant jumps) function aggregate(aggregator::Direct, u, p, t, end_time, constant_jumps, - ma_jumps, save_positions, rng; kwargs...) + ma_jumps, save_positions, rng; kwargs...) # handle constant jumps using tuples rates, affects! = get_jump_info_tuples(constant_jumps) build_jump_aggregation(DirectJumpAggregation, u, p, t, end_time, ma_jumps, - rates, affects!, save_positions, rng; kwargs...) + rates, affects!, save_positions, rng; kwargs...) end # creating the JumpAggregation structure (function wrapper-based constant jumps) function aggregate(aggregator::DirectFW, u, p, t, end_time, constant_jumps, - ma_jumps, save_positions, rng; kwargs...) + ma_jumps, save_positions, rng; kwargs...) # handle constant jumps using function wrappers rates, affects! = get_jump_info_fwrappers(u, p, t, constant_jumps) build_jump_aggregation(DirectJumpAggregation, u, p, t, end_time, ma_jumps, - rates, affects!, save_positions, rng; kwargs...) + rates, affects!, save_positions, rng; kwargs...) end # set up a new simulation and calculate the first jump / jump time @@ -53,7 +53,7 @@ end # execute one jump, changing the system state @inline function execute_jumps!(p::DirectJumpAggregation, integrator, u, params, t, - affects!) + affects!) update_state!(p, integrator, u, affects!) nothing end @@ -70,7 +70,7 @@ end # tuple-based constant jumps function time_to_next_jump(p::DirectJumpAggregation{T, S, F1}, u, params, - t) where {T, S, F1 <: Tuple} + t) where {T, S, F1 <: Tuple} prev_rate = zero(t) new_rate = zero(t) cur_rates = p.cur_rates @@ -112,7 +112,7 @@ end # function wrapper-based constant jumps function time_to_next_jump(p::DirectJumpAggregation{T, S, F1}, u, params, - t) where {T, S, F1 <: AbstractArray} + t) where {T, S, F1 <: AbstractArray} prev_rate = zero(t) new_rate = zero(t) cur_rates = p.cur_rates diff --git a/src/aggregators/directcr.jl b/src/aggregators/directcr.jl index b9aea31c..30d44216 100644 --- a/src/aggregators/directcr.jl +++ b/src/aggregators/directcr.jl @@ -11,7 +11,7 @@ by S. Mauch and M. Stalzer, ACM Trans. Comp. Biol. and Bioinf., 8, No. 1, 27-35 const MINJUMPRATE = 2.0^exponent(1e-12) mutable struct DirectCRJumpAggregation{T, S, F1, F2, RNG, DEPGR, U <: PriorityTable, - W <: Function} <: + W <: Function} <: AbstractSSAJumpAggregator{T, S, F1, F2, RNG} next_jump::Int prev_jump::Int @@ -32,11 +32,11 @@ mutable struct DirectCRJumpAggregation{T, S, F1, F2, RNG, DEPGR, U <: PriorityTa end function DirectCRJumpAggregation(nj::Int, njt::T, et::T, crs::Vector{T}, sr::T, - maj::S, rs::F1, affs!::F2, sps::Tuple{Bool, Bool}, - rng::RNG; num_specs, dep_graph = nothing, - minrate = convert(T, MINJUMPRATE), - maxrate = convert(T, Inf), - kwargs...) where {T, S, F1, F2, RNG} + maj::S, rs::F1, affs!::F2, sps::Tuple{Bool, Bool}, + rng::RNG; num_specs, dep_graph = nothing, + minrate = convert(T, MINJUMPRATE), + maxrate = convert(T, Inf), + kwargs...) where {T, S, F1, F2, RNG} # a dependency graph is needed and must be provided if there are constant rate jumps if dep_graph === nothing @@ -64,24 +64,24 @@ function DirectCRJumpAggregation(nj::Int, njt::T, et::T, crs::Vector{T}, sr::T, affecttype = F2 <: Tuple ? F2 : Any DirectCRJumpAggregation{T, S, F1, affecttype, RNG, typeof(dg), - typeof(rt), typeof(ratetogroup)}(nj, nj, njt, et, crs, sr, maj, - rs, affs!, sps, rng, dg, - minrate, maxrate, rt, - ratetogroup) + typeof(rt), typeof(ratetogroup)}(nj, nj, njt, et, crs, sr, maj, + rs, affs!, sps, rng, dg, + minrate, maxrate, rt, + ratetogroup) end ############################# Required Functions ############################## # creating the JumpAggregation structure (function wrapper-based constant jumps) function aggregate(aggregator::DirectCR, u, p, t, end_time, constant_jumps, - ma_jumps, save_positions, rng; kwargs...) + ma_jumps, save_positions, rng; kwargs...) # handle constant jumps using function wrappers rates, affects! = get_jump_info_fwrappers(u, p, t, constant_jumps) build_jump_aggregation(DirectCRJumpAggregation, u, p, t, end_time, ma_jumps, - rates, affects!, save_positions, rng; num_specs = length(u), - kwargs...) + rates, affects!, save_positions, rng; num_specs = length(u), + kwargs...) end # set up a new simulation and calculate the first jump / jump time diff --git a/src/aggregators/frm.jl b/src/aggregators/frm.jl index a4802e17..94baed37 100644 --- a/src/aggregators/frm.jl +++ b/src/aggregators/frm.jl @@ -13,35 +13,37 @@ mutable struct FRMJumpAggregation{T, S, F1, F2, RNG} <: rng::RNG end function FRMJumpAggregation(nj::Int, njt::T, et::T, crs::Vector{T}, sr::T, maj::S, rs::F1, - affs!::F2, sps::Tuple{Bool, Bool}, rng::RNG; - kwargs...) where {T, S, F1, F2, RNG} + affs!::F2, sps::Tuple{Bool, Bool}, rng::RNG; + kwargs...) where {T, S, F1, F2, RNG} affecttype = F2 <: Tuple ? F2 : Any FRMJumpAggregation{T, S, F1, affecttype, RNG}(nj, nj, njt, et, crs, sr, maj, rs, - affs!, sps, rng) + affs!, sps, rng) end ############################# Required Functions ############################# # creating the JumpAggregation structure (tuple-based constant jumps) function aggregate(aggregator::FRM, u, p, t, end_time, constant_jumps, - ma_jumps, save_positions, rng; kwargs...) + ma_jumps, save_positions, rng; kwargs...) # handle constant jumps using tuples rates, affects! = get_jump_info_tuples(constant_jumps) - build_jump_aggregation(FRMJumpAggregation, u, p, t, end_time, ma_jumps, rates, affects!, - save_positions, rng; kwargs...) + build_jump_aggregation( + FRMJumpAggregation, u, p, t, end_time, ma_jumps, rates, affects!, + save_positions, rng; kwargs...) end # creating the JumpAggregation structure (function wrapper-based constant jumps) function aggregate(aggregator::FRMFW, u, p, t, end_time, constant_jumps, - ma_jumps, save_positions, rng; kwargs...) + ma_jumps, save_positions, rng; kwargs...) # handle constant jumps using function wrappers rates, affects! = get_jump_info_fwrappers(u, p, t, constant_jumps) - build_jump_aggregation(FRMJumpAggregation, u, p, t, end_time, ma_jumps, rates, affects!, - save_positions, rng; kwargs...) + build_jump_aggregation( + FRMJumpAggregation, u, p, t, end_time, ma_jumps, rates, affects!, + save_positions, rng; kwargs...) end # set up a new simulation and calculate the first jump / jump time @@ -94,7 +96,7 @@ end # tuple-based constant jumps function next_constant_rate_jump(p::FRMJumpAggregation{T, S, F1, F2, RNG}, u, params, - t) where {T, S, F1 <: Tuple, F2 <: Tuple, RNG} + t) where {T, S, F1 <: Tuple, F2 <: Tuple, RNG} ttnj = typemax(typeof(t)) nextrx = zero(Int) if !isempty(p.rates) @@ -113,7 +115,7 @@ end # function wrapper-based constant jumps function next_constant_rate_jump(p::FRMJumpAggregation{T, S, F1}, u, params, - t) where {T, S, F1 <: AbstractArray} + t) where {T, S, F1 <: AbstractArray} ttnj = typemax(typeof(t)) nextrx = zero(Int) if !isempty(p.rates) diff --git a/src/aggregators/nrm.jl b/src/aggregators/nrm.jl index 8e086173..9f152f2f 100644 --- a/src/aggregators/nrm.jl +++ b/src/aggregators/nrm.jl @@ -19,9 +19,9 @@ mutable struct NRMJumpAggregation{T, S, F1, F2, RNG, DEPGR, PQ} <: end function NRMJumpAggregation(nj::Int, njt::T, et::T, crs::Vector{T}, sr::T, - maj::S, rs::F1, affs!::F2, sps::Tuple{Bool, Bool}, - rng::RNG; num_specs, dep_graph = nothing, - kwargs...) where {T, S, F1, F2, RNG} + maj::S, rs::F1, affs!::F2, sps::Tuple{Bool, Bool}, + rng::RNG; num_specs, dep_graph = nothing, + kwargs...) where {T, S, F1, F2, RNG} # a dependency graph is needed and must be provided if there are constant rate jumps if dep_graph === nothing @@ -41,22 +41,22 @@ function NRMJumpAggregation(nj::Int, njt::T, et::T, crs::Vector{T}, sr::T, affecttype = F2 <: Tuple ? F2 : Any NRMJumpAggregation{T, S, F1, affecttype, RNG, typeof(dg), typeof(pq)}(nj, nj, njt, et, - crs, sr, maj, - rs, affs!, sps, - rng, dg, pq) + crs, sr, maj, + rs, affs!, sps, + rng, dg, pq) end +############################# Required Functions ############################## # creating the JumpAggregation structure (function wrapper-based constant jumps) function aggregate(aggregator::NRM, u, p, t, end_time, constant_jumps, - ma_jumps, save_positions, rng; kwargs...) + ma_jumps, save_positions, rng; kwargs...) # handle constant jumps using function wrappers rates, affects! = get_jump_info_fwrappers(u, p, t, constant_jumps) build_jump_aggregation(NRMJumpAggregation, u, p, t, end_time, ma_jumps, - rates, affects!, save_positions, rng; num_specs = length(u), - kwargs...) + rates, affects!, save_positions, rng; num_specs = length(u), + kwargs...) end # set up a new simulation and calculate the first jump / jump time @@ -97,7 +97,7 @@ function update_dependent_rates!(p::NRMJumpAggregation, u, params, t) # update the jump rate @inbounds cur_rates[rx] = calculate_jump_rate(ma_jumps, num_majumps, rates, u, - params, t, rx) + params, t, rx) # calculate new jump times for dependent jumps if rx != p.next_jump && oldrate > zero(oldrate) diff --git a/src/aggregators/prioritytable.jl b/src/aggregators/prioritytable.jl index 9c394d37..0534dbf3 100644 --- a/src/aggregators/prioritytable.jl +++ b/src/aggregators/prioritytable.jl @@ -101,7 +101,7 @@ Setup table from a vector of priorities. The id of a priority is its position within this vector. """ function PriorityTable(priortogid::Function, priorities::AbstractVector, minpriority, - maxpriority) + maxpriority) numgroups = priortogid(maxpriority) numgroups -= one(typeof(numgroups)) pidtype = typeof(numgroups) @@ -348,7 +348,7 @@ end # DEFAULT NUMBINS: 20 * √length(times) # DEFAULT BINWIDTH: 16 / sum(propensities) function PriorityTimeTable(times::AbstractVector, mintime, timestep; binwidthconst = 16, - numbinsconst = 20) + numbinsconst = 20) binwidth = binwidthconst * timestep numbins = floor(Int64, numbinsconst * sqrt(length(times))) maxtime = mintime + numbins * binwidth @@ -365,7 +365,7 @@ function PriorityTimeTable(times::AbstractVector, mintime, timestep; binwidthcon end ptt = PriorityTimeTable(groups, pidtogroup, times, ttgdata, zero(pidtype), - zero(pidtype), maxtime, binwidthconst, numbinsconst) + zero(pidtype), maxtime, binwidthconst, numbinsconst) # Insert priority ids into the groups for (pid, time) in enumerate(times) if time > maxtime diff --git a/src/aggregators/rdirect.jl b/src/aggregators/rdirect.jl index 70534c17..e6949090 100644 --- a/src/aggregators/rdirect.jl +++ b/src/aggregators/rdirect.jl @@ -22,10 +22,10 @@ mutable struct RDirectJumpAggregation{T, S, F1, F2, RNG, DEPGR} <: end function RDirectJumpAggregation(nj::Int, njt::T, et::T, crs::Vector{T}, sr::T, maj::S, - rs::F1, affs!::F2, sps::Tuple{Bool, Bool}, rng::RNG; - num_specs, counter_threshold = length(crs), - dep_graph = nothing, - kwargs...) where {T, S, F1, F2, RNG} + rs::F1, affs!::F2, sps::Tuple{Bool, Bool}, rng::RNG; + num_specs, counter_threshold = length(crs), + dep_graph = nothing, + kwargs...) where {T, S, F1, F2, RNG} # a dependency graph is needed and must be provided if there are constant rate jumps if dep_graph === nothing if (get_num_majumps(maj) == 0) || !isempty(rs) @@ -43,24 +43,24 @@ function RDirectJumpAggregation(nj::Int, njt::T, et::T, crs::Vector{T}, sr::T, m max_rate = maximum(crs) affecttype = F2 <: Tuple ? F2 : Any return RDirectJumpAggregation{T, S, F1, affecttype, RNG, typeof(dg)}(nj, nj, njt, et, - crs, sr, maj, rs, - affs!, sps, rng, - dg, max_rate, 0, - counter_threshold) + crs, sr, maj, rs, + affs!, sps, rng, + dg, max_rate, 0, + counter_threshold) end ############################# Required Functions ############################# # creating the JumpAggregation structure (tuple-based constant jumps) function aggregate(aggregator::RDirect, u, p, t, end_time, constant_jumps, - ma_jumps, save_positions, rng; kwargs...) + ma_jumps, save_positions, rng; kwargs...) # handle constant jumps using function wrappers rates, affects! = get_jump_info_fwrappers(u, p, t, constant_jumps) build_jump_aggregation(RDirectJumpAggregation, u, p, t, end_time, ma_jumps, - rates, affects!, save_positions, rng; num_specs = length(u), - kwargs...) + rates, affects!, save_positions, rng; num_specs = length(u), + kwargs...) end # set up a new simulation and calculate the first jump / jump time @@ -116,8 +116,9 @@ function update_dependent_rates!(p::RDirectJumpAggregation, u, params, t) num_majumps = get_num_majumps(ma_jumps) max_rate_increased = false @inbounds for rx in dep_rxs - @inbounds new_rate = calculate_jump_rate(ma_jumps, num_majumps, rates, u, params, t, - rx) + @inbounds new_rate = calculate_jump_rate( + ma_jumps, num_majumps, rates, u, params, t, + rx) sum_rate += new_rate - cur_rates[rx] if new_rate > p.max_rate p.max_rate = new_rate diff --git a/src/aggregators/rssa.jl b/src/aggregators/rssa.jl index 8f6ed907..ee3ce19c 100644 --- a/src/aggregators/rssa.jl +++ b/src/aggregators/rssa.jl @@ -26,10 +26,10 @@ mutable struct RSSAJumpAggregation{T, S, F1, F2, RNG, VJMAP, JVMAP, BD, U} <: end function RSSAJumpAggregation(nj::Int, njt::T, et::T, crs::Vector{T}, sr::T, - maj::S, rs::F1, affs!::F2, sps::Tuple{Bool, Bool}, - rng::RNG; u::U, vartojumps_map = nothing, - jumptovars_map = nothing, - bracket_data = nothing, kwargs...) where {T, S, F1, F2, RNG, U} + maj::S, rs::F1, affs!::F2, sps::Tuple{Bool, Bool}, + rng::RNG; u::U, vartojumps_map = nothing, + jumptovars_map = nothing, + bracket_data = nothing, kwargs...) where {T, S, F1, F2, RNG, U} # a dependency graph is needed and must be provided if there are constant rate jumps if vartojumps_map === nothing if (get_num_majumps(maj) == 0) || !isempty(rs) @@ -64,24 +64,24 @@ function RSSAJumpAggregation(nj::Int, njt::T, et::T, crs::Vector{T}, sr::T, affecttype = F2 <: Tuple ? F2 : Any RSSAJumpAggregation{T, S, F1, affecttype, RNG, typeof(vtoj_map), - typeof(jtov_map), typeof(bd), U}(nj, nj, njt, et, crl_bnds, - crh_bnds, sr, maj, rs, affs!, sps, - rng, vtoj_map, jtov_map, bd, ulow, - uhigh) + typeof(jtov_map), typeof(bd), U}(nj, nj, njt, et, crl_bnds, + crh_bnds, sr, maj, rs, affs!, sps, + rng, vtoj_map, jtov_map, bd, ulow, + uhigh) end ############################# Required Functions ############################## # creating the JumpAggregation structure (function wrapper-based constant jumps) function aggregate(aggregator::RSSA, u, p, t, end_time, constant_jumps, - ma_jumps, save_positions, rng; kwargs...) + ma_jumps, save_positions, rng; kwargs...) # handle constant jumps using function wrappers rates, affects! = get_jump_info_fwrappers(u, p, t, constant_jumps) build_jump_aggregation(RSSAJumpAggregation, u, p, t, end_time, ma_jumps, - rates, affects!, save_positions, rng; u = u, - kwargs...) + rates, affects!, save_positions, rng; u = u, + kwargs...) end # set up a new simulation and calculate the first jump / jump time @@ -120,7 +120,7 @@ function generate_jumps!(p::RSSAJumpAggregation, integrator, u, params, t) end rerl += randexp(rng) @inbounds while rejectrx(ma_jumps, num_majumps, rates, cur_rate_high, - cur_rate_low, rng, u, jidx, params, t) + cur_rate_low, rng, u, jidx, params, t) # sample candidate reaction r = rand(rng) * sum_rate jidx = linear_search(cur_rate_high, r) diff --git a/src/aggregators/rssacr.jl b/src/aggregators/rssacr.jl index dc7f42a6..2f65e2f6 100644 --- a/src/aggregators/rssacr.jl +++ b/src/aggregators/rssacr.jl @@ -5,7 +5,7 @@ Composition-Rejection with Rejection sampling method (RSSA-CR) const MINJUMPRATE = 2.0^exponent(1e-12) mutable struct RSSACRJumpAggregation{F, S, F1, F2, RNG, U, VJMAP, JVMAP, BD, - P <: PriorityTable, W <: Function} <: + P <: PriorityTable, W <: Function} <: AbstractSSAJumpAggregator{F, S, F1, F2, RNG} next_jump::Int prev_jump::Int @@ -31,11 +31,11 @@ mutable struct RSSACRJumpAggregation{F, S, F1, F2, RNG, U, VJMAP, JVMAP, BD, end function RSSACRJumpAggregation(nj::Int, njt::F, et::F, crs::Vector{F}, sum_rate::F, maj::S, - rs::F1, affs!::F2, sps::Tuple{Bool, Bool}, rng::RNG; u::U, - vartojumps_map = nothing, jumptovars_map = nothing, - bracket_data = nothing, minrate = convert(F, MINJUMPRATE), - maxrate = convert(F, Inf), - kwargs...) where {F, S, F1, F2, RNG, U} + rs::F1, affs!::F2, sps::Tuple{Bool, Bool}, rng::RNG; u::U, + vartojumps_map = nothing, jumptovars_map = nothing, + bracket_data = nothing, minrate = convert(F, MINJUMPRATE), + maxrate = convert(F, Inf), + kwargs...) where {F, S, F1, F2, RNG, U} # a dependency graph is needed and must be provided if there are constant rate jumps if vartojumps_map === nothing if (get_num_majumps(maj) == 0) || !isempty(rs) @@ -81,24 +81,24 @@ function RSSACRJumpAggregation(nj::Int, njt::F, et::F, crs::Vector{F}, sum_rate: affecttype = F2 <: Tuple ? F2 : Any RSSACRJumpAggregation{typeof(njt), S, F1, affecttype, RNG, U, typeof(vtoj_map), - typeof(jtov_map), typeof(bd), typeof(rt), - typeof(ratetogroup)}(nj, nj, njt, et, crl_bnds, crh_bnds, - sum_rate, maj, rs, affs!, sps, rng, vtoj_map, - jtov_map, bd, ulow, uhigh, minrate, maxrate, - rt, ratetogroup) + typeof(jtov_map), typeof(bd), typeof(rt), + typeof(ratetogroup)}(nj, nj, njt, et, crl_bnds, crh_bnds, + sum_rate, maj, rs, affs!, sps, rng, vtoj_map, + jtov_map, bd, ulow, uhigh, minrate, maxrate, + rt, ratetogroup) end ############################# Required Functions ############################## # creating the JumpAggregation structure (function wrapper-based constant jumps) function aggregate(aggregator::RSSACR, u, p, t, end_time, constant_jumps, - ma_jumps, save_positions, rng; kwargs...) + ma_jumps, save_positions, rng; kwargs...) # handle constant jumps using function wrappers rates, affects! = get_jump_info_fwrappers(u, p, t, constant_jumps) build_jump_aggregation(RSSACRJumpAggregation, u, p, t, end_time, ma_jumps, - rates, affects!, save_positions, rng; u = u, kwargs...) + rates, affects!, save_positions, rng; u = u, kwargs...) end # set up a new simulation and calculate the first jump / jump time @@ -145,7 +145,7 @@ function generate_jumps!(p::RSSACRJumpAggregation, integrator, u, params, t) end rerl += randexp(rng) while rejectrx(ma_jumps, num_majumps, rates, cur_rate_high, cur_rate_low, rng, u, jidx, - params, t) + params, t) # sample candidate reaction jidx = sample(rt, cur_rate_high, rng) rerl += randexp(rng) @@ -162,7 +162,7 @@ end update bracketing for species that depend on the just executed jump """ @inline function update_dependent_rates!(p::RSSACRJumpAggregation, u::AbstractVector, - params, t) + params, t) # update bracketing intervals @unpack ulow, uhigh = p crhigh = p.cur_rate_high diff --git a/src/aggregators/sortingdirect.jl b/src/aggregators/sortingdirect.jl index 078ac3c9..f9048f03 100644 --- a/src/aggregators/sortingdirect.jl +++ b/src/aggregators/sortingdirect.jl @@ -21,9 +21,9 @@ mutable struct SortingDirectJumpAggregation{T, S, F1, F2, RNG, DEPGR} <: end function SortingDirectJumpAggregation(nj::Int, njt::T, et::T, crs::Vector{T}, sr::T, - maj::S, rs::F1, affs!::F2, sps::Tuple{Bool, Bool}, - rng::RNG; num_specs, dep_graph = nothing, - kwargs...) where {T, S, F1, F2, RNG} + maj::S, rs::F1, affs!::F2, sps::Tuple{Bool, Bool}, + rng::RNG; num_specs, dep_graph = nothing, + kwargs...) where {T, S, F1, F2, RNG} # a dependency graph is needed and must be provided if there are constant rate jumps if dep_graph === nothing @@ -43,24 +43,24 @@ function SortingDirectJumpAggregation(nj::Int, njt::T, et::T, crs::Vector{T}, sr jtoidx = collect(1:length(crs)) affecttype = F2 <: Tuple ? F2 : Any SortingDirectJumpAggregation{T, S, F1, affecttype, RNG, typeof(dg)}(nj, nj, njt, et, - crs, sr, maj, rs, - affs!, sps, rng, - dg, jtoidx, - zero(Int)) + crs, sr, maj, rs, + affs!, sps, rng, + dg, jtoidx, + zero(Int)) end ############################# Required Functions ############################## # creating the JumpAggregation structure (function wrapper-based constant jumps) function aggregate(aggregator::SortingDirect, u, p, t, end_time, constant_jumps, - ma_jumps, save_positions, rng; kwargs...) + ma_jumps, save_positions, rng; kwargs...) # handle constant jumps using function wrappers rates, affects! = get_jump_info_fwrappers(u, p, t, constant_jumps) build_jump_aggregation(SortingDirectJumpAggregation, u, p, t, end_time, ma_jumps, - rates, affects!, save_positions, rng; num_specs = length(u), - kwargs...) + rates, affects!, save_positions, rng; num_specs = length(u), + kwargs...) end # set up a new simulation and calculate the first jump / jump time diff --git a/src/aggregators/ssajump.jl b/src/aggregators/ssajump.jl index 6dd51909..3f6df60a 100644 --- a/src/aggregators/ssajump.jl +++ b/src/aggregators/ssajump.jl @@ -46,7 +46,7 @@ end end @inline function concretize_affects!(p::AbstractSSAJumpAggregator, - ::I) where {I <: DiffEqBase.DEIntegrator} + ::I) where {I <: DiffEqBase.DEIntegrator} if (p.affects! isa Vector) && !(p.affects! isa Vector{FunctionWrappers.FunctionWrapper{Nothing, Tuple{I}}}) AffectWrapper = FunctionWrappers.FunctionWrapper{Nothing, Tuple{I}} @@ -56,8 +56,8 @@ end end @inline function concretize_affects!(p::AbstractSSAJumpAggregator{T, S, F1, F2}, - ::I) where {T, S, F1, F2 <: Tuple, - I <: DiffEqBase.DEIntegrator} + ::I) where {T, S, F1, F2 <: Tuple, + I <: DiffEqBase.DEIntegrator} nothing end @@ -89,8 +89,8 @@ function (p::AbstractSSAJumpAggregator)(integrator::I) where {I <: DiffEqBase.DE end function (p::AbstractSSAJumpAggregator{ - T, S, F1, F2})(integrator::DiffEqBase.DEIntegrator) where - {T, S, F1, F2 <: Union{Tuple, Nothing}} + T, S, F1, F2})(integrator::DiffEqBase.DEIntegrator) where + {T, S, F1, F2 <: Union{Tuple, Nothing}} execute_jumps!(p, integrator, integrator.u, integrator.p, integrator.t, p.affects!) generate_jumps!(p, integrator, integrator.u, integrator.p, integrator.t) register_next_jump_time!(integrator, p, integrator.t) @@ -118,14 +118,14 @@ end Helper routine for setting up standard fields of SSA jump aggregations. """ function build_jump_aggregation(jump_agg_type, u, p, t, end_time, ma_jumps, rates, - affects!, save_positions, rng; kwargs...) + affects!, save_positions, rng; kwargs...) # mass action jumps majumps = ma_jumps if majumps === nothing majumps = MassActionJump(Vector{typeof(t)}(), - Vector{Vector{Pair{Int, eltype(u)}}}(), - Vector{Vector{Pair{Int, eltype(u)}}}()) + Vector{Vector{Pair{Int, eltype(u)}}}(), + Vector{Vector{Pair{Int, eltype(u)}}}()) end # current jump rates, allows mass action rates and constant jumps @@ -135,7 +135,7 @@ function build_jump_aggregation(jump_agg_type, u, p, t, end_time, ma_jumps, rate next_jump = 0 next_jump_time = typemax(typeof(t)) jump_agg_type(next_jump, next_jump_time, end_time, cur_rates, sum_rate, - majumps, rates, affects!, save_positions, rng; kwargs...) + majumps, rates, affects!, save_positions, rng; kwargs...) end """ @@ -197,7 +197,7 @@ function update_dependent_rates!(p::AbstractSSAJumpAggregator, u, params, t) @inbounds for rx in dep_rxs sum_rate -= cur_rates[rx] @inbounds cur_rates[rx] = calculate_jump_rate(p.ma_jumps, num_majumps, p.rates, u, - params, t, rx) + params, t, rx) sum_rate += cur_rates[rx] end @@ -230,7 +230,7 @@ Execute `p.next_jump`. end @generated function update_state!(p::AbstractSSAJumpAggregator, integrator, u, - affects!::T) where {T <: Tuple} + affects!::T) where {T <: Tuple} quote @unpack ma_jumps, next_jump = p num_ma_rates = get_num_majumps(ma_jumps) @@ -297,8 +297,9 @@ end Perform rejection sampling test (used in RSSA methods). """ -@inline function rejectrx(ma_jumps, num_majumps, rates, cur_rate_high, cur_rate_low, rng, u, - jidx, params, t) +@inline function rejectrx( + ma_jumps, num_majumps, rates, cur_rate_high, cur_rate_low, rng, u, + jidx, params, t) # rejection test @inbounds r2 = rand(rng) * cur_rate_high[jidx] @inbounds crlow = cur_rate_low[jidx] diff --git a/src/coupling.jl b/src/coupling.jl index eae21794..c3417b2e 100644 --- a/src/coupling.jl +++ b/src/coupling.jl @@ -4,11 +4,11 @@ methods for stochastically modeled population processes. IMA J Numer Anal 2015; 35 (4): 1757-1778. doi: 10.1093/imanum/dru044 """ function SplitCoupledJumpProblem(prob::DiffEqBase.AbstractJumpProblem, - prob_control::DiffEqBase.AbstractJumpProblem, - aggregator::AbstractAggregatorAlgorithm, - coupling_map::Vector{Tuple{Int, Int}}; kwargs...) + prob_control::DiffEqBase.AbstractJumpProblem, + aggregator::AbstractAggregatorAlgorithm, + coupling_map::Vector{Tuple{Int, Int}}; kwargs...) JumpProblem(cat_problems(prob.prob, prob_control.prob), aggregator, - build_split_jumps(prob, prob_control, coupling_map)...; kwargs...) + build_split_jumps(prob, prob_control, coupling_map)...; kwargs...) end # make new problem by joining initial_data @@ -18,7 +18,7 @@ function cat_problems(prob::DiscreteProblem, prob_control::DiscreteProblem) end function cat_problems(prob::DiffEqBase.AbstractODEProblem, - prob_control::DiffEqBase.AbstractODEProblem) + prob_control::DiffEqBase.AbstractODEProblem) l = length(prob.u0) # add l_c = length(prob_control.u0) _f = SciMLBase.unwrapped_f(prob.f) @@ -50,7 +50,7 @@ function cat_problems(prob::DiscreteProblem, prob_control::DiffEqBase.AbstractOD end function cat_problems(prob::DiffEqBase.AbstractSDEProblem, - prob_control::DiffEqBase.AbstractSDEProblem) + prob_control::DiffEqBase.AbstractSDEProblem) l = length(prob.u0) new_f = function (du, u, p, t) prob.f(@view(du[1:l]), u.u, p, t) @@ -65,7 +65,7 @@ function cat_problems(prob::DiffEqBase.AbstractSDEProblem, end function cat_problems(prob::DiffEqBase.AbstractSDEProblem, - prob_control::DiffEqBase.AbstractODEProblem) + prob_control::DiffEqBase.AbstractODEProblem) l = length(prob.u0) _f = SciMLBase.unwrapped_f(prob.f) @@ -111,14 +111,14 @@ function cat_problems(prob_control::DiscreteProblem, prob::DiffEqBase.AbstractSD cat_problems(prob, prob_control) end function cat_problems(prob_control::DiffEqBase.AbstractODEProblem, - prob::DiffEqBase.AbstractSDEProblem) + prob::DiffEqBase.AbstractSDEProblem) cat_problems(prob, prob_control) end # this only depends on the jumps in prob, not prob.prob function build_split_jumps(prob::DiffEqBase.AbstractJumpProblem, - prob_control::DiffEqBase.AbstractJumpProblem, - coupling_map::Vector{Tuple{Int, Int}}) + prob_control::DiffEqBase.AbstractJumpProblem, + coupling_map::Vector{Tuple{Int, Int}}) num_jumps = length(prob.discrete_jump_aggregation.rates) num_jumps_control = length(prob_control.discrete_jump_aggregation.rates) jumps = [] diff --git a/src/extended_jump_array.jl b/src/extended_jump_array.jl index 20db5b7d..fe43692e 100644 --- a/src/extended_jump_array.jl +++ b/src/extended_jump_array.jl @@ -96,7 +96,7 @@ end # Ignore axes function Base.similar(A::ExtendedJumpArray, ::Type{S}, - axes::Tuple{Base.OneTo{Int}}) where {S} + axes::Tuple{Base.OneTo{Int}}) where {S} ExtendedJumpArray(similar(A.u, S), similar(A.jump_u, S)) end @@ -126,49 +126,49 @@ plot_indices(A::ExtendedJumpArray) = eachindex(A.u) # The jump array styles stores two sub-styles in the type, # one for the `u` array and one for the `jump_u` array struct ExtendedJumpArrayStyle{UStyle <: Broadcast.BroadcastStyle, - JumpUStyle <: Broadcast.BroadcastStyle} <: + JumpUStyle <: Broadcast.BroadcastStyle} <: Broadcast.BroadcastStyle end # Init style based on type of u/jump_u function ExtendedJumpArrayStyle(::US, ::JumpUS) where {US, JumpUS} ExtendedJumpArrayStyle{US, JumpUS}() end function Base.BroadcastStyle(::Type{ExtendedJumpArray{ - T3, T1, UType, JumpUType}}) where {T3, - T1, - UType, - JumpUType - } + T3, T1, UType, JumpUType}}) where {T3, + T1, + UType, + JumpUType +} ExtendedJumpArrayStyle(Base.BroadcastStyle(UType), Base.BroadcastStyle(JumpUType)) end # Combine with other styles by combining individually with u/jump_u styles function Base.BroadcastStyle(::ExtendedJumpArrayStyle{UStyle, JumpUStyle}, - ::Style) where {UStyle, JumpUStyle, - Style <: Base.Broadcast.BroadcastStyle} + ::Style) where {UStyle, JumpUStyle, + Style <: Base.Broadcast.BroadcastStyle} ExtendedJumpArrayStyle(Broadcast.result_style(UStyle(), Style()), - Broadcast.result_style(JumpUStyle(), Style())) + Broadcast.result_style(JumpUStyle(), Style())) end # Decay back to the DefaultArrayStyle for higher-order default styles, to support adding to raw vectors as needed function Base.BroadcastStyle(::ExtendedJumpArrayStyle{UStyle, JumpUStyle}, - ::Broadcast.DefaultArrayStyle{0}) where {UStyle, JumpUStyle} + ::Broadcast.DefaultArrayStyle{0}) where {UStyle, JumpUStyle} ExtendedJumpArrayStyle(UStyle(), JumpUStyle()) end function Base.BroadcastStyle(::ExtendedJumpArrayStyle{UStyle, JumpUStyle}, - ::Broadcast.DefaultArrayStyle{N}) where {N, UStyle, JumpUStyle} + ::Broadcast.DefaultArrayStyle{N}) where {N, UStyle, JumpUStyle} Broadcast.DefaultArrayStyle{N}() end function Base.Broadcast.BroadcastStyle(::S, - ::Base.Broadcast.Unknown) where { - UStyle, JumpUStyle, - S <: - JumpProcesses.ExtendedJumpArrayStyle{ - UStyle, - JumpUStyle - } - } + ::Base.Broadcast.Unknown) where { + UStyle, JumpUStyle, + S <: + JumpProcesses.ExtendedJumpArrayStyle{ + UStyle, + JumpUStyle + } +} return throw(ArgumentError("Cannot broadcast JumpProcesses.ExtendedJumpArray with" * " something of type Base.Broadcast.Unknown.")) end @@ -185,7 +185,7 @@ find_eja(a::ExtendedJumpArray, rest) = a find_eja(::Any, rest) = find_eja(rest) function Base.similar(bc::Broadcast.Broadcasted{ExtendedJumpArrayStyle{US, JumpUS}}, - ::Type{ElType}) where {US, JumpUS, ElType} + ::Type{ElType}) where {US, JumpUS, ElType} A = find_eja(bc) ExtendedJumpArray(similar(A.u, ElType), similar(A.jump_u, ElType)) end @@ -195,11 +195,11 @@ end Broadcast.Broadcasted{Style}(bc.f, repack_args(i, bc.args)) end @inline function repack(bc::Broadcast.Broadcasted{ExtendedJumpArrayStyle{US, JumpUS}}, - i::Val{:u}) where {US, JumpUS} + i::Val{:u}) where {US, JumpUS} Broadcast.Broadcasted{US}(bc.f, repack_args(i, bc.args)) end @inline function repack(bc::Broadcast.Broadcasted{ExtendedJumpArrayStyle{US, JumpUS}}, - i::Val{:jump_u}) where {US, JumpUS} + i::Val{:jump_u}) where {US, JumpUS} Broadcast.Broadcasted{JumpUS}(bc.f, repack_args(i, bc.args)) end @@ -219,10 +219,10 @@ end end @inline function Base.copyto!(dest::ExtendedJumpArray, - bc::Broadcast.Broadcasted{ExtendedJumpArrayStyle{US, JumpUS}}) where { - US, - JumpUS - } + bc::Broadcast.Broadcasted{ExtendedJumpArrayStyle{US, JumpUS}}) where { + US, + JumpUS +} copyto!(dest.u, repack(bc, Val(:u))) copyto!(dest.jump_u, repack(bc, Val(:jump_u))) dest diff --git a/src/jumps.jl b/src/jumps.jl index a990c738..910f2175 100644 --- a/src/jumps.jl +++ b/src/jumps.jl @@ -172,12 +172,12 @@ function VariableRateJump(rate, affect!; lrate = nothing, urate = nothing, ``` """ function VariableRateJump(rate, affect!; - lrate = nothing, urate = nothing, - rateinterval = nothing, rootfind = true, - idxs = nothing, - save_positions = (false, true), - interp_points = 10, - abstol = 1e-12, reltol = 0) + lrate = nothing, urate = nothing, + rateinterval = nothing, rootfind = true, + idxs = nothing, + save_positions = (false, true), + interp_points = 10, + abstol = 1e-12, reltol = 0) if !(urate !== nothing && rateinterval !== nothing) && !(urate === nothing && rateinterval === nothing) error("`urate` and `rateinterval` must both be `nothing`, or must both be defined.") @@ -189,7 +189,7 @@ function VariableRateJump(rate, affect!; end VariableRateJump(rate, affect!, lrate, urate, rateinterval, idxs, rootfind, - interp_points, save_positions, abstol, reltol) + interp_points, save_positions, abstol, reltol) end """ @@ -337,8 +337,8 @@ struct MassActionJump{T, S, U, V} <: AbstractMassActionJump param_mapper::V function MassActionJump{T, S, U, V}(rates::T, rs_in::S, ns::U, pmapper::V, - scale_rates::Bool, useiszero::Bool, - nocopy::Bool) where {T <: AbstractVector, S, U, V} + scale_rates::Bool, useiszero::Bool, + nocopy::Bool) where {T <: AbstractVector, S, U, V} sr = nocopy ? rates : copy(rates) rs = nocopy ? rs_in : copy(rs_in) for i in eachindex(rs) @@ -353,12 +353,12 @@ struct MassActionJump{T, S, U, V} <: AbstractMassActionJump new(sr, rs, ns, pmapper) end function MassActionJump{Nothing, Vector{S}, - Vector{U}, V}(::Nothing, rs_in::Vector{S}, - ns::Vector{U}, pmapper::V, - scale_rates::Bool, - useiszero::Bool, - nocopy::Bool) where {S <: AbstractVector, - U <: AbstractVector, V} + Vector{U}, V}(::Nothing, rs_in::Vector{S}, + ns::Vector{U}, pmapper::V, + scale_rates::Bool, + useiszero::Bool, + nocopy::Bool) where {S <: AbstractVector, + U <: AbstractVector, V} rs = nocopy ? rs_in : copy(rs_in) for i in eachindex(rs) if useiszero && (length(rs[i]) == 1) && iszero(rs[i][1][1]) @@ -368,8 +368,8 @@ struct MassActionJump{T, S, U, V} <: AbstractMassActionJump new(nothing, rs, ns, pmapper) end function MassActionJump{T, S, U, V}(rate::T, rs_in::S, ns::U, pmapper::V, - scale_rates::Bool, useiszero::Bool, - nocopy::Bool) where {T <: Number, S, U, V} + scale_rates::Bool, useiszero::Bool, + nocopy::Bool) where {T <: Number, S, U, V} rs = rs_in if useiszero && (length(rs) == 1) && iszero(rs[1][1]) rs = typeof(rs)() @@ -378,8 +378,8 @@ struct MassActionJump{T, S, U, V} <: AbstractMassActionJump new(sr, rs, ns, pmapper) end function MassActionJump{Nothing, S, U, V}(::Nothing, rs_in::S, ns::U, pmapper::V, - scale_rates::Bool, useiszero::Bool, - nocopy::Bool) where {S, U, V} + scale_rates::Bool, useiszero::Bool, + nocopy::Bool) where {S, U, V} rs = rs_in if useiszero && (length(rs) == 1) && iszero(rs[1][1]) rs = typeof(rs)() @@ -388,23 +388,23 @@ struct MassActionJump{T, S, U, V} <: AbstractMassActionJump end end function MassActionJump(usr::T, rs::S, ns::U, pmapper::V; scale_rates = true, - useiszero = true, nocopy = false) where {T, S, U, V} + useiszero = true, nocopy = false) where {T, S, U, V} MassActionJump{T, S, U, V}(usr, rs, ns, pmapper, scale_rates, useiszero, nocopy) end function MassActionJump(usr::T, rs, ns; scale_rates = true, useiszero = true, - nocopy = false) where {T <: AbstractVector} + nocopy = false) where {T <: AbstractVector} MassActionJump(usr, rs, ns, nothing; scale_rates = scale_rates, useiszero = useiszero, - nocopy = nocopy) + nocopy = nocopy) end function MassActionJump(usr::T, rs, ns; scale_rates = true, useiszero = true, - nocopy = false) where {T <: Number} + nocopy = false) where {T <: Number} MassActionJump(usr, rs, ns, nothing; scale_rates = scale_rates, useiszero = useiszero, - nocopy = nocopy) + nocopy = nocopy) end # with parameter indices or mapping, multiple jump case function MassActionJump(rs, ns; param_idxs = nothing, param_mapper = nothing, - scale_rates = true, useiszero = true, nocopy = false) + scale_rates = true, useiszero = true, nocopy = false) if param_mapper === nothing (param_idxs === nothing) && error("If no parameter indices are given via param_idxs, an explicit parameter mapping must be passed in via param_mapper.") @@ -416,7 +416,7 @@ function MassActionJump(rs, ns; param_idxs = nothing, param_mapper = nothing, end MassActionJump(nothing, nocopy ? rs : copy(rs), ns, pmapper; scale_rates = scale_rates, - useiszero = useiszero, nocopy = true) + useiszero = useiszero, nocopy = true) end using_params(maj::MassActionJump{T, S, U, Nothing}) where {T, S, U} = false @@ -442,8 +442,8 @@ end # update a maj with parameter vectors function (ratemap::MassActionJumpParamMapper{U})(maj::MassActionJump, newparams; - scale_rates, - kwargs...) where {U <: AbstractArray} + scale_rates, + kwargs...) where {U <: AbstractArray} for i in 1:get_num_majumps(maj) maj.scaled_rates[i] = newparams[ratemap.param_idxs[i]] end @@ -456,18 +456,18 @@ function to_collection(ratemap::MassActionJumpParamMapper{Int}) end function Base.merge!(pmap1::MassActionJumpParamMapper{U}, - pmap2::MassActionJumpParamMapper{U}) where {U <: AbstractVector} + pmap2::MassActionJumpParamMapper{U}) where {U <: AbstractVector} append!(pmap1.param_idxs, pmap2.param_idxs) end function Base.merge!(pmap1::MassActionJumpParamMapper{U}, - pmap2::MassActionJumpParamMapper{V}) where {U <: AbstractVector, - V <: Int} + pmap2::MassActionJumpParamMapper{V}) where {U <: AbstractVector, + V <: Int} push!(pmap1.param_idxs, pmap2.param_idxs) end function Base.merge(pmap1::MassActionJumpParamMapper{Int}, - pmap2::MassActionJumpParamMapper{Int}) + pmap2::MassActionJumpParamMapper{Int}) MassActionJumpParamMapper([pmap1.param_idxs, pmap2.param_idxs]) end @@ -542,7 +542,7 @@ JumpSet(jump::VariableRateJump) = JumpSet((jump,), (), nothing, nothing) JumpSet(jump::RegularJump) = JumpSet((), (), jump, nothing) JumpSet(jump::AbstractMassActionJump) = JumpSet((), (), nothing, jump) function JumpSet(; variable_jumps = (), constant_jumps = (), - regular_jumps = nothing, massaction_jumps = nothing) + regular_jumps = nothing, massaction_jumps = nothing) JumpSet(variable_jumps, constant_jumps, regular_jumps, massaction_jumps) end JumpSet(jb::Nothing) = JumpSet() @@ -559,7 +559,7 @@ function JumpSet(vjs, cjs, rj, majv::Vector{T}) where {T <: MassActionJump} end maj = setup_majump_to_merge(majv[1].scaled_rates, majv[1].reactant_stoch, - majv[1].net_stoch, majv[1].param_mapper) + majv[1].net_stoch, majv[1].param_mapper) for i in 2:length(majv) massaction_jump_combine(maj, majv[i]) end @@ -605,9 +605,9 @@ end end @inline function split_jumps(vj, cj, rj, maj, j::JumpSet, args...) split_jumps((vj..., j.variable_jumps...), - (cj..., j.constant_jumps...), - regular_jump_combine(rj, j.regular_jump), - massaction_jump_combine(maj, j.massaction_jump), args...) + (cj..., j.constant_jumps...), + regular_jump_combine(rj, j.regular_jump), + massaction_jump_combine(maj, j.massaction_jump), args...) end regular_jump_combine(rj1::RegularJump, rj2::Nothing) = rj1 @@ -620,42 +620,42 @@ end # functionality to merge two mass action jumps together function check_majump_type(maj::MassActionJump{S, T, U, V}) where {S <: Number, T, U, V} setup_majump_to_merge(maj.scaled_rates, maj.reactant_stoch, maj.net_stoch, - maj.param_mapper) + maj.param_mapper) end function check_majump_type(maj::MassActionJump{Nothing, T, U, V}) where {T, U, V} setup_majump_to_merge(maj.scaled_rates, maj.reactant_stoch, maj.net_stoch, - maj.param_mapper) + maj.param_mapper) end # if given containers of rates and stoichiometry directly create a jump function setup_majump_to_merge(sr::T, rs::AbstractVector{S}, ns::AbstractVector{U}, - pmapper) where {T <: AbstractVector, S <: AbstractArray, - U <: AbstractArray} + pmapper) where {T <: AbstractVector, S <: AbstractArray, + U <: AbstractArray} MassActionJump(sr, rs, ns, pmapper; scale_rates = false) end # if just given the data for one jump (and not in a container) wrap in a vector function setup_majump_to_merge(sr::S, rs::T, ns::U, - pmapper) where {S <: Number, T <: AbstractArray, - U <: AbstractArray} + pmapper) where {S <: Number, T <: AbstractArray, + U <: AbstractArray} MassActionJump([sr], [rs], [ns], - (pmapper === nothing) ? pmapper : to_collection(pmapper); - scale_rates = false) + (pmapper === nothing) ? pmapper : to_collection(pmapper); + scale_rates = false) end # if no rate field setup yet function setup_majump_to_merge(::Nothing, rs::T, ns::U, - pmapper) where {T <: AbstractArray, U <: AbstractArray} + pmapper) where {T <: AbstractArray, U <: AbstractArray} MassActionJump(nothing, [rs], [ns], - (pmapper === nothing) ? pmapper : to_collection(pmapper); - scale_rates = false) + (pmapper === nothing) ? pmapper : to_collection(pmapper); + scale_rates = false) end # when given a collection of reactions to add to maj function majump_merge!(maj::MassActionJump{U, <:AbstractVector{V}, <:AbstractVector{W}, X}, - sr::U, rs::AbstractVector{V}, ns::AbstractVector{W}, - param_mapper) where {U <: Union{AbstractVector, Nothing}, - V <: AbstractVector, W <: AbstractVector, X} + sr::U, rs::AbstractVector{V}, ns::AbstractVector{W}, + param_mapper) where {U <: Union{AbstractVector, Nothing}, + V <: AbstractVector, W <: AbstractVector, X} (U <: AbstractVector) && append!(maj.scaled_rates, sr) append!(maj.reactant_stoch, rs) append!(maj.net_stoch, ns) @@ -670,11 +670,11 @@ end # when given a single jump's worth of data to add to maj function majump_merge!(maj::MassActionJump{U, V, W, X}, sr::T, rs::S1, ns::S2, - param_mapper) where {T <: Union{Number, Nothing}, - S1 <: AbstractArray, S2 <: AbstractArray, - U <: Union{AbstractVector{T}, Nothing}, - V <: AbstractVector{S1}, - W <: AbstractVector{S2}, X} + param_mapper) where {T <: Union{Number, Nothing}, + S1 <: AbstractArray, S2 <: AbstractArray, + U <: Union{AbstractVector{T}, Nothing}, + V <: AbstractVector{S1}, + W <: AbstractVector{S2}, X} (T <: Number) && push!(maj.scaled_rates, sr) push!(maj.reactant_stoch, rs) push!(maj.net_stoch, ns) @@ -691,18 +691,18 @@ end # when maj only stores a single jump's worth of data (and not in a collection) # create a new jump with the merged data stored in vectors function majump_merge!(maj::MassActionJump{T, S, U, V}, sr::T, rs::S, ns::U, - param_mapper::V) where {T <: Union{Number, Nothing}, - S <: AbstractArray{<:Pair}, - U <: AbstractArray{<:Pair}, V} + param_mapper::V) where {T <: Union{Number, Nothing}, + S <: AbstractArray{<:Pair}, + U <: AbstractArray{<:Pair}, V} rates = (T <: Nothing) ? nothing : [maj.scaled_rates, sr] if maj.param_mapper === nothing (param_mapper === nothing) || error("Error, trying to merge a MassActionJump with a parameter mapping to one without a parameter mapping.") return MassActionJump(rates, [maj.reactant_stoch, rs], [maj.net_stoch, ns], - param_mapper; scale_rates = false) + param_mapper; scale_rates = false) else return MassActionJump(rates, [maj.reactant_stoch, rs], [maj.net_stoch, ns], - merge(maj.param_mapper, param_mapper); scale_rates = false) + merge(maj.param_mapper, param_mapper); scale_rates = false) end end @@ -711,7 +711,7 @@ massaction_jump_combine(maj1::Nothing, maj2::MassActionJump) = maj2 massaction_jump_combine(maj1::Nothing, maj2::Nothing) = maj1 function massaction_jump_combine(maj1::MassActionJump, maj2::MassActionJump) majump_merge!(maj1, maj2.scaled_rates, maj2.reactant_stoch, maj2.net_stoch, - maj2.param_mapper) + maj2.param_mapper) end ##### helper methods for unpacking rates and affects! from constant jumps ##### @@ -729,7 +729,7 @@ end function get_jump_info_fwrappers(u, p, t, constant_jumps) RateWrapper = FunctionWrappers.FunctionWrapper{typeof(t), - Tuple{typeof(u), typeof(p), typeof(t)}} + Tuple{typeof(u), typeof(p), typeof(t)}} if (constant_jumps !== nothing) && !isempty(constant_jumps) rates = [RateWrapper(c.rate) for c in constant_jumps] diff --git a/src/massaction_rates.jl b/src/massaction_rates.jl index 17273368..d44bfff6 100644 --- a/src/massaction_rates.jl +++ b/src/massaction_rates.jl @@ -4,8 +4,8 @@ ############################################################################### @inline function evalrxrate(speciesvec::AbstractVector{T}, rxidx::S, - majump::MassActionJump{U, V, W, X})::R where - {T, S, R, U <: AbstractVector{R}, V, W, X} + majump::MassActionJump{U, V, W, X})::R where + {T, S, R, U <: AbstractVector{R}, V, W, X} val = one(T) @inbounds for specstoch in majump.reactant_stoch[rxidx] specpop = speciesvec[specstoch[1]] @@ -20,7 +20,7 @@ end @inline function executerx!(speciesvec::AbstractVector{T}, rxidx::S, - majump::M) where {T, S, M <: AbstractMassActionJump} + majump::M) where {T, S, M <: AbstractMassActionJump} @inbounds net_stoch = majump.net_stoch[rxidx] @inbounds for specstoch in net_stoch speciesvec[specstoch[1]] += specstoch[2] @@ -29,11 +29,11 @@ end end @inline function executerx(speciesvec::SVector{T}, rxidx::S, - majump::M) where {T, S, M <: AbstractMassActionJump} + majump::M) where {T, S, M <: AbstractMassActionJump} @inbounds net_stoch = majump.net_stoch[rxidx] @inbounds for specstoch in net_stoch speciesvec = setindex(speciesvec, speciesvec[specstoch[1]] + specstoch[2], - specstoch[1]) + specstoch[1]) end speciesvec @@ -45,8 +45,8 @@ end end function scalerates!(unscaled_rates::AbstractVector{U}, - stochmat::AbstractVector{V}) where {U, S, T, W <: Pair{S, T}, - V <: AbstractVector{W}} + stochmat::AbstractVector{V}) where {U, S, T, W <: Pair{S, T}, + V <: AbstractVector{W}} @inbounds for i in eachindex(unscaled_rates) coef = one(T) @inbounds for specstoch in stochmat[i] @@ -58,8 +58,8 @@ function scalerates!(unscaled_rates::AbstractVector{U}, end function scalerates!(unscaled_rates::AbstractMatrix{U}, - stochmat::AbstractVector{V}) where {U, S, T, W <: Pair{S, T}, - V <: AbstractVector{W}} + stochmat::AbstractVector{V}) where {U, S, T, W <: Pair{S, T}, + V <: AbstractVector{W}} @inbounds for i in size(unscaled_rates, 1) coef = one(T) @inbounds for specstoch in stochmat[i] @@ -71,7 +71,7 @@ function scalerates!(unscaled_rates::AbstractMatrix{U}, end function scalerate(unscaled_rate::U, - stochmat::AbstractVector{Pair{S, T}}) where {U <: Number, S, T} + stochmat::AbstractVector{Pair{S, T}}) where {U <: Number, S, T} coef = one(T) @inbounds for specstoch in stochmat coef *= factorial(specstoch[2]) diff --git a/src/problem.jl b/src/problem.jl index 38d27b3e..e684789f 100644 --- a/src/problem.jl +++ b/src/problem.jl @@ -65,7 +65,7 @@ DifferentialEquations.jl [docs](https://docs.sciml.ai/JumpProcesses/stable/) for commonly asked questions. """ mutable struct JumpProblem{iip, P, A, C, J <: Union{Nothing, AbstractJumpAggregator}, J2, - J3, J4, R, K} <: DiffEqBase.AbstractJumpProblem{P, J} + J3, J4, R, K} <: DiffEqBase.AbstractJumpProblem{P, J} """The type of problem to couple the jumps to. For a pure jump process use `DiscreteProblem`, to couple to ODEs, `ODEProblem`, etc.""" prob::P """The aggregator algorithm that determines the next jump times and types for `ConstantRateJump`s and `MassActionJump`s. Examples include `Direct`.""" @@ -86,7 +86,7 @@ mutable struct JumpProblem{iip, P, A, C, J <: Union{Nothing, AbstractJumpAggrega kwargs::K end function JumpProblem(p::P, a::A, dj::J, jc::C, vj::J2, rj::J3, mj::J4, - rng::R, kwargs::K) where {P, A, J, C, J2, J3, J4, R, K} + rng::R, kwargs::K) where {P, A, J, C, J2, J3, J4, R, K} iip = isinplace_jump(p, rj) JumpProblem{iip, P, A, C, J, J2, J3, J4, R, K}(p, a, dj, jc, vj, rj, mj, rng, kwargs) end @@ -151,8 +151,8 @@ function DiffEqBase.remake(jprob::JumpProblem; kwargs...) end T(newprob, jprob.aggregator, jprob.discrete_jump_aggregation, jprob.jump_callback, - jprob.variable_jumps, jprob.regular_jump, jprob.massaction_jump, jprob.rng, - jprob.kwargs) + jprob.variable_jumps, jprob.regular_jump, jprob.massaction_jump, jprob.rng, + jprob.kwargs) end # when setindex! is used. @@ -193,28 +193,28 @@ function JumpProblem(prob, jumps::AbstractJump...; kwargs...) end function JumpProblem(prob, aggregator::AbstractAggregatorAlgorithm, - jumps::ConstantRateJump; kwargs...) + jumps::ConstantRateJump; kwargs...) JumpProblem(prob, aggregator, JumpSet(jumps); kwargs...) end function JumpProblem(prob, aggregator::AbstractAggregatorAlgorithm, - jumps::VariableRateJump; kwargs...) + jumps::VariableRateJump; kwargs...) JumpProblem(prob, aggregator, JumpSet(jumps); kwargs...) end function JumpProblem(prob, aggregator::AbstractAggregatorAlgorithm, jumps::RegularJump; - kwargs...) + kwargs...) JumpProblem(prob, aggregator, JumpSet(jumps); kwargs...) end function JumpProblem(prob, aggregator::AbstractAggregatorAlgorithm, - jumps::AbstractMassActionJump; kwargs...) + jumps::AbstractMassActionJump; kwargs...) JumpProblem(prob, aggregator, JumpSet(jumps); kwargs...) end function JumpProblem(prob, aggregator::AbstractAggregatorAlgorithm, jumps::AbstractJump...; - kwargs...) + kwargs...) JumpProblem(prob, aggregator, JumpSet(jumps...); kwargs...) end function JumpProblem(prob, jumps::JumpSet; vartojumps_map = nothing, - jumptovars_map = nothing, dep_graph = nothing, - spatial_system = nothing, hopping_constants = nothing, kwargs...) + jumptovars_map = nothing, dep_graph = nothing, + spatial_system = nothing, hopping_constants = nothing, kwargs...) ps = (; vartojumps_map, jumptovars_map, dep_graph, spatial_system, hopping_constants) aggtype = select_aggregator(jumps; ps...) return JumpProblem(prob, aggtype(), jumps; ps..., kwargs...) @@ -228,20 +228,20 @@ end make_kwarg(; kwargs...) = kwargs function JumpProblem(prob, aggregator::AbstractAggregatorAlgorithm, jumps::JumpSet; - save_positions = prob isa DiffEqBase.AbstractDiscreteProblem ? - (false, true) : (true, true), - rng = DEFAULT_RNG, scale_rates = true, useiszero = true, - spatial_system = nothing, hopping_constants = nothing, - callback = nothing, use_vrj_bounds = true, kwargs...) + save_positions = prob isa DiffEqBase.AbstractDiscreteProblem ? + (false, true) : (true, true), + rng = DEFAULT_RNG, scale_rates = true, useiszero = true, + spatial_system = nothing, hopping_constants = nothing, + callback = nothing, use_vrj_bounds = true, kwargs...) # initialize the MassActionJump rate constants with the user parameters if using_params(jumps.massaction_jump) rates = jumps.massaction_jump.param_mapper(prob.p) maj = MassActionJump(rates, jumps.massaction_jump.reactant_stoch, - jumps.massaction_jump.net_stoch, - jumps.massaction_jump.param_mapper; scale_rates = scale_rates, - useiszero = useiszero, - nocopy = true) + jumps.massaction_jump.net_stoch, + jumps.massaction_jump.param_mapper; scale_rates = scale_rates, + useiszero = useiszero, + nocopy = true) else maj = jumps.massaction_jump end @@ -279,7 +279,7 @@ function JumpProblem(prob, aggregator::AbstractAggregatorAlgorithm, jumps::JumpS constant_jump_callback = CallbackSet() else disc_agg = aggregate(aggregator, u, prob.p, t, end_time, jumps.constant_jumps, maj, - save_positions, rng; kwargs...) + save_positions, rng; kwargs...) constant_jump_callback = DiscreteCallback(disc_agg) end @@ -299,13 +299,13 @@ function JumpProblem(prob, aggregator::AbstractAggregatorAlgorithm, jumps::JumpS solkwargs = make_kwarg(; callback) JumpProblem{iip, typeof(new_prob), typeof(aggregator), - typeof(jump_cbs), typeof(disc_agg), - typeof(cont_agg), - typeof(jumps.regular_jump), - typeof(maj), typeof(rng), typeof(solkwargs)}(new_prob, aggregator, disc_agg, - jump_cbs, cont_agg, - jumps.regular_jump, maj, rng, - solkwargs) + typeof(jump_cbs), typeof(disc_agg), + typeof(cont_agg), + typeof(jumps.regular_jump), + typeof(maj), typeof(rng), typeof(solkwargs)}(new_prob, aggregator, disc_agg, + jump_cbs, cont_agg, + jumps.regular_jump, maj, rng, + solkwargs) end # extends prob.u0 to an ExtendedJumpArray with Njumps integrated intensity values, @@ -342,7 +342,7 @@ function extend_problem(prob::DiffEqBase.AbstractODEProblem, jumps; rng = DEFAUL u0 = extend_u0(prob, length(jumps), rng) f = ODEFunction{isinplace(prob)}(jump_f; sys = prob.f.sys, - observed = prob.f.observed) + observed = prob.f.observed) remake(prob; f, u0) end @@ -378,7 +378,7 @@ function extend_problem(prob::DiffEqBase.AbstractSDEProblem, jumps; rng = DEFAUL u0 = extend_u0(prob, length(jumps), rng) f = SDEFunction{isinplace(prob)}(jump_f, jump_g; sys = prob.f.sys, - observed = prob.f.observed) + observed = prob.f.observed) remake(prob; f, g = jump_g, u0) end @@ -404,7 +404,7 @@ function extend_problem(prob::DiffEqBase.AbstractDDEProblem, jumps; rng = DEFAUL u0 = extend_u0(prob, length(jumps), rng) f = DDEFunction{isinplace(prob)}(jump_f; sys = prob.f.sys, - observed = prob.f.observed) + observed = prob.f.observed) remake(prob; f, u0) end @@ -431,7 +431,7 @@ function extend_problem(prob::DiffEqBase.AbstractDAEProblem, jumps; rng = DEFAUL u0 = extend_u0(prob, length(jumps), rng) f = DAEFunction{isinplace(prob)}(jump_f, sys = prob.f.sys, - observed = prob.f.observed) + observed = prob.f.observed) remake(prob; f, u0) end @@ -445,12 +445,12 @@ function build_variable_callback(cb, idx, jump, jumps...; rng = DEFAULT_RNG) integrator.u.jump_u[idx] = -randexp(rng, typeof(integrator.t)) end new_cb = ContinuousCallback(condition, affect!; - idxs = jump.idxs, - rootfind = jump.rootfind, - interp_points = jump.interp_points, - save_positions = jump.save_positions, - abstol = jump.abstol, - reltol = jump.reltol) + idxs = jump.idxs, + rootfind = jump.rootfind, + interp_points = jump.interp_points, + save_positions = jump.save_positions, + abstol = jump.abstol, + reltol = jump.reltol) build_variable_callback(CallbackSet(cb, new_cb), idx, jumps...; rng = DEFAULT_RNG) end @@ -464,19 +464,19 @@ function build_variable_callback(cb, idx, jump; rng = DEFAULT_RNG) integrator.u.jump_u[idx] = -randexp(rng, typeof(integrator.t)) end new_cb = ContinuousCallback(condition, affect!; - idxs = jump.idxs, - rootfind = jump.rootfind, - interp_points = jump.interp_points, - save_positions = jump.save_positions, - abstol = jump.abstol, - reltol = jump.reltol) + idxs = jump.idxs, + rootfind = jump.rootfind, + interp_points = jump.interp_points, + save_positions = jump.save_positions, + abstol = jump.abstol, + reltol = jump.reltol) CallbackSet(cb, new_cb) end aggregator(jp::JumpProblem{iip, P, A, C, J}) where {iip, P, A, C, J} = A @inline function extend_tstops!(tstops, - jp::JumpProblem{P, A, C, J, J2}) where {P, A, C, J, J2} + jp::JumpProblem{P, A, C, J, J2}) where {P, A, C, J, J2} !(jp.jump_callback.discrete_callbacks isa Tuple{}) && push!(tstops, jp.jump_callback.discrete_callbacks[1].condition.next_jump_time) end @@ -498,18 +498,18 @@ num_constant_rate_jumps(aggregator::AbstractSSAJumpAggregator) = length(aggregat function Base.summary(io::IO, prob::JumpProblem) type_color, no_color = SciMLBase.get_colorizers(io) print(io, - type_color, nameof(typeof(prob)), - no_color, " with problem ", - type_color, nameof(typeof(prob.prob)), - no_color, " with aggregator ", - type_color, typeof(prob.aggregator)) + type_color, nameof(typeof(prob)), + no_color, " with problem ", + type_color, nameof(typeof(prob.prob)), + no_color, " with aggregator ", + type_color, typeof(prob.aggregator)) end function Base.show(io::IO, mime::MIME"text/plain", A::JumpProblem) summary(io, A) println(io) println(io, "Number of jumps with discrete aggregation: ", - A.discrete_jump_aggregation === nothing ? 0 : - num_constant_rate_jumps(A.discrete_jump_aggregation)) + A.discrete_jump_aggregation === nothing ? 0 : + num_constant_rate_jumps(A.discrete_jump_aggregation)) println(io, "Number of jumps with continuous aggregation: ", length(A.variable_jumps)) nmajs = (A.massaction_jump !== nothing) ? get_num_majumps(A.massaction_jump) : 0 println(io, "Number of mass action jumps: ", nmajs) diff --git a/src/simple_regular_solve.jl b/src/simple_regular_solve.jl index 2ebbb15d..5cc58f49 100644 --- a/src/simple_regular_solve.jl +++ b/src/simple_regular_solve.jl @@ -1,8 +1,8 @@ struct SimpleTauLeaping <: DiffEqBase.DEAlgorithm end function DiffEqBase.solve(jump_prob::JumpProblem, alg::SimpleTauLeaping; - seed = nothing, - dt = error("dt is required for SimpleTauLeaping.")) + seed = nothing, + dt = error("dt is required for SimpleTauLeaping.")) # boilerplate from SimpleTauLeaping method @assert isempty(jump_prob.jump_callback.continuous_callbacks) # still needs to be a regular jump @@ -46,8 +46,8 @@ function DiffEqBase.solve(jump_prob::JumpProblem, alg::SimpleTauLeaping; end sol = DiffEqBase.build_solution(prob, alg, t, u, - calculate_error = false, - interp = DiffEqBase.ConstantInterpolation(t, u)) + calculate_error = false, + interp = DiffEqBase.ConstantInterpolation(t, u)) end export SimpleTauLeaping diff --git a/src/solve.jl b/src/solve.jl index cb646ada..ee51ce98 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -1,7 +1,7 @@ function DiffEqBase.__solve(jump_prob::DiffEqBase.AbstractJumpProblem{P}, - alg::DiffEqBase.DEAlgorithm, timeseries = [], ts = [], ks = [], - recompile::Type{Val{recompile_flag}} = Val{true}; - kwargs...) where {P, recompile_flag} + alg::DiffEqBase.DEAlgorithm, timeseries = [], ts = [], ks = [], + recompile::Type{Val{recompile_flag}} = Val{true}; + kwargs...) where {P, recompile_flag} integrator = init(jump_prob, alg, timeseries, ts, ks, recompile; kwargs...) solve!(integrator) integrator.sol @@ -10,7 +10,7 @@ end # if passed a JumpProblem over a DiscreteProblem, and no aggregator is selected use # SSAStepper function DiffEqBase.__solve(jump_prob::DiffEqBase.AbstractJumpProblem{P}; - kwargs...) where {P <: DiscreteProblem} + kwargs...) where {P <: DiscreteProblem} DiffEqBase.__solve(jump_prob, SSAStepper(); kwargs...) end @@ -19,11 +19,11 @@ function DiffEqBase.__solve(jump_prob::DiffEqBase.AbstractJumpProblem; kwargs... end function DiffEqBase.__init(_jump_prob::DiffEqBase.AbstractJumpProblem{P}, - alg::DiffEqBase.DEAlgorithm, timeseries = [], ts = [], ks = [], - recompile::Type{Val{recompile_flag}} = Val{true}; - callback = nothing, seed = nothing, - alias_jump = Threads.threadid() == 1, - kwargs...) where {P, recompile_flag} + alg::DiffEqBase.DEAlgorithm, timeseries = [], ts = [], ks = [], + recompile::Type{Val{recompile_flag}} = Val{true}; + callback = nothing, seed = nothing, + alias_jump = Threads.threadid() == 1, + kwargs...) where {P, recompile_flag} if alias_jump jump_prob = _jump_prob reset_jump_problem!(jump_prob, seed) @@ -35,13 +35,13 @@ function DiffEqBase.__init(_jump_prob::DiffEqBase.AbstractJumpProblem{P}, if jump_prob.prob isa DiffEqBase.AbstractDDEProblem # callback comes after jump consistent with SSAStepper integrator = init(jump_prob.prob, alg, timeseries, ts, ks; - callback = CallbackSet(jump_prob.jump_callback, callback), - kwargs...) + callback = CallbackSet(jump_prob.jump_callback, callback), + kwargs...) else # callback comes after jump consistent with SSAStepper integrator = init(jump_prob.prob, alg, timeseries, ts, ks, recompile; - callback = CallbackSet(jump_prob.jump_callback, callback), - kwargs...) + callback = CallbackSet(jump_prob.jump_callback, callback), + kwargs...) end end @@ -50,7 +50,7 @@ function resetted_jump_problem(_jump_prob, seed) if !isempty(jump_prob.jump_callback.discrete_callbacks) if seed === nothing Random.seed!(jump_prob.jump_callback.discrete_callbacks[1].condition.rng, - rand(UInt64)) + rand(UInt64)) else Random.seed!(jump_prob.jump_callback.discrete_callbacks[1].condition.rng, seed) end diff --git a/src/spatial/bracketing.jl b/src/spatial/bracketing.jl index 7ef147e0..b8ac2095 100644 --- a/src/spatial/bracketing.jl +++ b/src/spatial/bracketing.jl @@ -35,7 +35,7 @@ end function total_site_rate(rx_rates::LowHigh, hop_rates::LowHigh, site) return LowHigh(total_site_rate(rx_rates.low, hop_rates.low, site), - total_site_rate(rx_rates.high, hop_rates.high, site)) + total_site_rate(rx_rates.high, hop_rates.high, site)) end function update_rx_rates!(rx_rates::LowHigh, rxs, u_low_high, integrator, site) diff --git a/src/spatial/directcrdirect.jl b/src/spatial/directcrdirect.jl index 4ad7b0a9..bc4f54cb 100644 --- a/src/spatial/directcrdirect.jl +++ b/src/spatial/directcrdirect.jl @@ -6,8 +6,8 @@ const MINJUMPRATE = 2.0^exponent(1e-12) #NOTE state vector u is a matrix. u[i,j] is species i, site j #NOTE hopping_constants is a matrix. hopping_constants[i,j] is species i, site j mutable struct DirectCRDirectJumpAggregation{T, S, F1, F2, RNG, J, RX, HOP, DEPGR, - VJMAP, JVMAP, SS, U <: PriorityTable, - W <: Function} <: + VJMAP, JVMAP, SS, U <: PriorityTable, + W <: Function} <: AbstractSSAJumpAggregator{T, S, F1, F2, RNG} next_jump::SpatialJump{J} #some structure to identify the next event: reaction or hop prev_jump::SpatialJump{J} #some structure to identify the previous event: reaction or hop @@ -30,12 +30,12 @@ mutable struct DirectCRDirectJumpAggregation{T, S, F1, F2, RNG, J, RX, HOP, DEPG end function DirectCRDirectJumpAggregation(nj::SpatialJump{J}, njt::T, et::T, rx_rates::RX, - hop_rates::HOP, site_rates::Vector{T}, - sps::Tuple{Bool, Bool}, rng::RNG, spatial_system::SS; - num_specs, minrate = convert(T, MINJUMPRATE), - vartojumps_map = nothing, jumptovars_map = nothing, - dep_graph = nothing, - kwargs...) where {J, T, RX, HOP, RNG, SS} + hop_rates::HOP, site_rates::Vector{T}, + sps::Tuple{Bool, Bool}, rng::RNG, spatial_system::SS; + num_specs, minrate = convert(T, MINJUMPRATE), + vartojumps_map = nothing, jumptovars_map = nothing, + dep_graph = nothing, + kwargs...) where {J, T, RX, HOP, RNG, SS} # a dependency graph is needed if dep_graph === nothing @@ -70,26 +70,26 @@ function DirectCRDirectJumpAggregation(nj::SpatialJump{J}, njt::T, et::T, rx_rat rt = PriorityTable(ratetogroup, zeros(T, 1), minrate, 2 * minrate) DirectCRDirectJumpAggregation{T, Nothing, Nothing, Nothing, RNG, J, RX, HOP, - typeof(dg), typeof(vtoj_map), - typeof(jtov_map), SS, typeof(rt), - typeof(ratetogroup)}(nj, nj, njt, et, rx_rates, hop_rates, - site_rates, nothing, nothing, sps, - rng, dg, vtoj_map, - jtov_map, spatial_system, num_specs, - rt, ratetogroup) + typeof(dg), typeof(vtoj_map), + typeof(jtov_map), SS, typeof(rt), + typeof(ratetogroup)}(nj, nj, njt, et, rx_rates, hop_rates, + site_rates, nothing, nothing, sps, + rng, dg, vtoj_map, + jtov_map, spatial_system, num_specs, + rt, ratetogroup) end ############################# Required Functions ############################## # creating the JumpAggregation structure (function wrapper-based constant jumps) function aggregate(aggregator::DirectCRDirect, starting_state, p, t, end_time, - constant_jumps, ma_jumps, save_positions, rng; hopping_constants, - spatial_system, kwargs...) + constant_jumps, ma_jumps, save_positions, rng; hopping_constants, + spatial_system, kwargs...) num_species = size(starting_state, 1) majumps = ma_jumps if majumps === nothing majumps = MassActionJump(Vector{typeof(end_time)}(), - Vector{Vector{Pair{Int, Int}}}(), - Vector{Vector{Pair{Int, Int}}}()) + Vector{Vector{Pair{Int, Int}}}(), + Vector{Vector{Pair{Int, Int}}}()) end next_jump = SpatialJump{Int}(typemax(Int), typemax(Int), typemax(Int)) #a placeholder @@ -99,8 +99,8 @@ function aggregate(aggregator::DirectCRDirect, starting_state, p, t, end_time, site_rates = zeros(typeof(end_time), num_sites(spatial_system)) DirectCRDirectJumpAggregation(next_jump, next_jump_time, end_time, rx_rates, hop_rates, - site_rates, save_positions, rng, spatial_system; - num_specs = num_species, kwargs...) + site_rates, save_positions, rng, spatial_system; + num_specs = num_species, kwargs...) end # set up a new simulation and calculate the first jump / jump time @@ -122,7 +122,7 @@ end # execute one jump, changing the system state function execute_jumps!(p::DirectCRDirectJumpAggregation, integrator, u, params, t, - affects!) + affects!) # execute jump update_state!(p, integrator) @@ -138,7 +138,7 @@ end reset all structs, reevaluate all rates, repopulate the priority table """ function fill_rates_and_get_times!(aggregation::DirectCRDirectJumpAggregation, integrator, - t) + t) @unpack spatial_system, rx_rates, hop_rates, site_rates, rt = aggregation u = integrator.u @@ -168,7 +168,7 @@ end recalculate jump rates for jumps that depend on the just executed jump (p.prev_jump) """ function update_dependent_rates_and_firing_times!(p::DirectCRDirectJumpAggregation, - integrator, t) + integrator, t) u = integrator.u site_rates = p.site_rates jump = p.prev_jump diff --git a/src/spatial/flatten.jl b/src/spatial/flatten.jl index 4a5e7a1b..390b0e32 100644 --- a/src/spatial/flatten.jl +++ b/src/spatial/flatten.jl @@ -4,13 +4,13 @@ using JumpProcesses, DiffEqBase, Graphs prob.u0 must be a Matrix with prob.u0[i,j] being the number of species i at site j """ function flatten(ma_jump, prob::DiscreteProblem, spatial_system, hopping_constants; - kwargs...) + kwargs...) tspan = prob.tspan u0 = prob.u0 if ma_jump === nothing ma_jump = MassActionJump(Vector{typeof(tspan[1])}(), - Vector{Vector{Pair{Int, eltype(u0)}}}(), - Vector{Vector{Pair{Int, eltype(u0)}}}()) + Vector{Vector{Pair{Int, eltype(u0)}}}(), + Vector{Vector{Pair{Int, eltype(u0)}}}()) end netstoch = ma_jump.net_stoch reactstoch = ma_jump.reactant_stoch @@ -24,25 +24,25 @@ function flatten(ma_jump, prob::DiscreteProblem, spatial_system, hopping_constan rx_rates = ma_jump.spatial_rates elseif isnothing(ma_jump.spatial_rates) rx_rates = reshape(repeat(ma_jump.uniform_rates, num_nodes), - length(ma_jump.uniform_rates), num_nodes) + length(ma_jump.uniform_rates), num_nodes) else @assert size(ma_jump.spatial_rates, 2) == num_nodes rx_rates = cat(dims = 1, - reshape(repeat(ma_jump.uniform_rates, num_nodes), - length(ma_jump.uniform_rates), num_nodes), - ma_jump.spatial_rates) + reshape(repeat(ma_jump.uniform_rates, num_nodes), + length(ma_jump.uniform_rates), num_nodes), + ma_jump.spatial_rates) end end flatten(netstoch, reactstoch, rx_rates, spatial_system, u0, tspan, hopping_constants; - scale_rates = false, kwargs...) + scale_rates = false, kwargs...) end """ if hopping_constants is a matrix, assume hopping_constants[i,j] is the hopping constant of species i from site j to any neighbor """ function flatten(netstoch::AbstractArray, reactstoch::AbstractArray, - rx_rates::AbstractArray, spatial_system, u0::Matrix{Int}, tspan, - hopping_constants::Matrix{F}; kwargs...) where {F <: Number} + rx_rates::AbstractArray, spatial_system, u0::Matrix{Int}, tspan, + hopping_constants::Matrix{F}; kwargs...) where {F <: Number} @assert size(hopping_constants) == size(u0) hop_constants = Matrix{Vector{F}}(undef, size(hopping_constants)) for ci in CartesianIndices(hop_constants) @@ -51,28 +51,28 @@ function flatten(netstoch::AbstractArray, reactstoch::AbstractArray, ones(outdegree(spatial_system, site)) end flatten(netstoch, reactstoch, rx_rates, spatial_system, u0, tspan, hop_constants; - kwargs...) + kwargs...) end """ if reaction rates is a vector, assume reaction rates are equal across sites """ function flatten(netstoch::AbstractArray, reactstoch::AbstractArray, rx_rates::Vector, - spatial_system, u0::Matrix{Int}, tspan, - hopping_constants::Matrix{Vector{F}}; kwargs...) where {F <: Number} + spatial_system, u0::Matrix{Int}, tspan, + hopping_constants::Matrix{Vector{F}}; kwargs...) where {F <: Number} num_nodes = num_sites(spatial_system) rates = reshape(repeat(rx_rates, num_nodes), length(rx_rates), num_nodes) flatten(netstoch, reactstoch, rates, spatial_system, u0, tspan, hopping_constants; - kwargs...) + kwargs...) end """ "flatten" the spatial jump problem. Return flattened DiscreteProblem and MassActionJump. """ function flatten(netstoch::Vector{R}, reactstoch::Vector{R}, rx_rates::Matrix{F}, - spatial_system, u0::Matrix{Int}, tspan, - hopping_constants::Matrix{Vector{F}}; scale_rates = true, - kwargs...) where {R, F <: Number} + spatial_system, u0::Matrix{Int}, tspan, + hopping_constants::Matrix{Vector{F}}; scale_rates = true, + kwargs...) where {R, F <: Number} num_species = size(u0, 1) num_nodes = num_sites(spatial_system) num_rxs = length(reactstoch) @@ -119,7 +119,7 @@ function flatten(netstoch::Vector{R}, reactstoch::Vector{R}, rx_rates::Matrix{F} # put everything together ma_jump = MassActionJump(total_rates, total_reactstoch, total_netstoch; nocopy = true, - scale_rates = scale_rates) + scale_rates = scale_rates) flattened_u0 = vec(u0) prob = DiscreteProblem(flattened_u0, tspan, total_rates) prob, ma_jump diff --git a/src/spatial/hop_rates.jl b/src/spatial/hop_rates.jl index 0b04e3e2..46cec2d5 100644 --- a/src/spatial/hop_rates.jl +++ b/src/spatial/hop_rates.jl @@ -14,7 +14,7 @@ function HopRates(hopping_constants::Vector{F}, spatial_system) where {F <: Numb HopRatesGraphDs(hopping_constants, num_sites(spatial_system)) end function HopRates(hopping_constants::Vector{F}, - grid::CartesianGridRej) where {F <: Number} + grid::CartesianGridRej) where {F <: Number} HopRatesGraphDs(hopping_constants, num_sites(grid)) end @@ -22,7 +22,7 @@ function HopRates(hopping_constants::Matrix{F}, spatial_system) where {F <: Numb HopRatesGraphDsi(hopping_constants) end function HopRates(hopping_constants::Matrix{F}, - grid::CartesianGridRej) where {F <: Number} + grid::CartesianGridRej) where {F <: Number} HopRatesGraphDsi(hopping_constants) end @@ -30,29 +30,29 @@ function HopRates(hopping_constants::Matrix{Vector{F}}, spatial_system) where {F HopRatesGraphDsij(hopping_constants) end function HopRates(hopping_constants::Matrix{Vector{F}}, - grid::CartesianGridRej) where {F <: Number} + grid::CartesianGridRej) where {F <: Number} HopRatesGridDsij(hopping_constants, grid) end function HopRates(p::Pair{SpecHop, SiteHop}, - spatial_system) where {F <: Number, SpecHop <: Vector{F}, - SiteHop <: Vector{Vector{F}}} + spatial_system) where {F <: Number, SpecHop <: Vector{F}, + SiteHop <: Vector{Vector{F}}} HopRatesGraphDsLij(p...) end function HopRates(p::Pair{SpecHop, SiteHop}, - grid::CartesianGridRej) where - {F <: Number, SpecHop <: Vector{F}, SiteHop <: Vector{Vector{F}}} + grid::CartesianGridRej) where + {F <: Number, SpecHop <: Vector{F}, SiteHop <: Vector{Vector{F}}} HopRatesGridDsLij(p..., grid) end function HopRates(p::Pair{SpecHop, SiteHop}, - spatial_system) where {F <: Number, SpecHop <: Matrix{F}, - SiteHop <: Vector{Vector{F}}} + spatial_system) where {F <: Number, SpecHop <: Matrix{F}, + SiteHop <: Vector{Vector{F}}} HopRatesGraphDsiLij(p...) end function HopRates(p::Pair{SpecHop, SiteHop}, - grid::CartesianGridRej) where - {SpecHop <: Matrix{F}, SiteHop <: Vector{Vector{F}}} where {F <: Number} + grid::CartesianGridRej) where + {SpecHop <: Matrix{F}, SiteHop <: Vector{Vector{F}}} where {F <: Number} HopRatesGridDsiLij(p..., grid) end @@ -62,7 +62,7 @@ end update rates of all specs in species at site """ function update_hop_rates!(hop_rates::AbstractHopRates, species::AbstractArray, u, site, - spatial_system) + spatial_system) @inbounds for spec in species update_hop_rate!(hop_rates, spec, u, site, spatial_system) end @@ -77,7 +77,7 @@ function update_hop_rate!(hop_rates::AbstractHopRates, species, u, site, spatial rates = hop_rates.rates @inbounds old_rate = rates[species, site] @inbounds rates[species, site] = evalhoprate(hop_rates, u, species, site, - spatial_system) + spatial_system) @inbounds hop_rates.sum_rates[site] += rates[species, site] - old_rate old_rate end @@ -107,7 +107,7 @@ sample species to hop from site """ function sample_species(hop_rates::AbstractHopRates, site, rng) @inbounds linear_search((@view hop_rates.rates[:, site]), - rand(rng) * total_site_hop_rate(hop_rates, site)) + rand(rng) * total_site_hop_rate(hop_rates, site)) end """ @@ -136,7 +136,7 @@ end function Base.show(io::IO, ::MIME"text/plain", hop_rates::HopRatesGraphDs) num_specs, num_sites = size(hop_rates.rates) println(io, - "HopRates with $num_specs species and $num_sites sites. \nHopping constants of form D_{s} where s is species.") + "HopRates with $num_specs species and $num_sites sites. \nHopping constants of form D_{s} where s is species.") end """ @@ -176,7 +176,7 @@ end function Base.show(io::IO, ::MIME"text/plain", hop_rates::HopRatesGraphDsi) num_specs, num_sites = size(hop_rates.rates) println(io, - "HopRates with $num_specs species and $num_sites sites. \nHopping constants of form D_{s,i} where s is species, and i is source.") + "HopRates with $num_specs species and $num_sites sites. \nHopping constants of form D_{s,i} where s is species, and i is source.") end """ @@ -216,7 +216,7 @@ end function Base.show(io::IO, ::MIME"text/plain", hop_rates::HopRatesGraphDsij) num_specs, num_sites = size(hop_rates.rates) println(io, - "HopRates with $num_specs species and $num_sites sites. \nHopping constants of form D_{s,i,j} where s is species, i is source and j is destination.") + "HopRates with $num_specs species and $num_sites sites. \nHopping constants of form D_{s,i,j} where s is species, i is source and j is destination.") end """ @@ -225,7 +225,7 @@ end initializes HopRates with zero rates """ function HopRatesGraphDsij(hopping_constants::Matrix{Vector{F}}; - do_cumsum = true) where {F <: Number} + do_cumsum = true) where {F <: Number} do_cumsum && (hopping_constants = map(cumsum, hopping_constants)) rates = zeros(F, size(hopping_constants)) sum_rates = zeros(F, size(rates, 2)) @@ -233,7 +233,7 @@ function HopRatesGraphDsij(hopping_constants::Matrix{Vector{F}}; end function sample_target_site(hop_rates::HopRatesGraphDsij, site, species, rng, - spatial_system) + spatial_system) @inbounds cum_hop_consts = hop_rates.hop_const_cumulative_sums[species, site] @inbounds n = searchsortedfirst(cum_hop_consts, rand(rng) * cum_hop_consts[end]) return nth_nbr(spatial_system, site, n) @@ -261,7 +261,7 @@ end function Base.show(io::IO, ::MIME"text/plain", hop_rates::HopRatesGridDsij) num_specs, num_sites = size(hop_rates.rates) println(io, - "HopRates with $num_specs species and $num_sites sites, optimized for CartesianGrid. \nHopping constants of form L_{s,i,j} where s is species, i is source and j is destination.") + "HopRates with $num_specs species and $num_sites sites, optimized for CartesianGrid. \nHopping constants of form L_{s,i,j} where s is species, i is source and j is destination.") end """ @@ -270,7 +270,7 @@ end initializes HopRates with zero rates """ function HopRatesGridDsij(hopping_constants::Array{F, 3}; - do_cumsum = true) where {F <: Number} + do_cumsum = true) where {F <: Number} do_cumsum && (hopping_constants = mapslices(cumsum, hopping_constants, dims = 1)) rates = zeros(F, size(hopping_constants)[2:3]) sum_rates = zeros(F, size(rates, 2)) @@ -279,7 +279,7 @@ end function HopRatesGridDsij(hopping_constants::Matrix{Vector{F}}, grid) where {F <: Number} new_hopping_constants = Array{F, 3}(undef, 2 * dimension(grid), - size(hopping_constants)...) + size(hopping_constants)...) for ci in CartesianIndices(hopping_constants) species, site = Tuple(ci) nb_constants = @view new_hopping_constants[:, species, site] @@ -317,9 +317,9 @@ end function Base.show(io::IO, ::MIME"text/plain", hop_rates::HopRatesGraphDsLij) num_specs, num_sites = length(hop_rates.species_hop_constants), - length(hop_rates.hop_const_cumulative_sums) + length(hop_rates.hop_const_cumulative_sums) println(io, - "HopRates with $num_specs species and $num_sites sites. \nHopping constants of form D_s * L_{i,j} where s is species, i is source and j is destination.") + "HopRates with $num_specs species and $num_sites sites. \nHopping constants of form D_s * L_{i,j} where s is species, i is source and j is destination.") end """ @@ -328,8 +328,8 @@ end initializes HopRates with zero rates """ function HopRatesGraphDsLij(species_hop_constants::Vector{F}, - site_hop_constants::Vector{Vector{F}}; - do_cumsum = true) where {F <: Number} + site_hop_constants::Vector{Vector{F}}; + do_cumsum = true) where {F <: Number} do_cumsum && (site_hop_constants = map(cumsum, site_hop_constants)) rates = zeros(F, length(species_hop_constants), length(site_hop_constants)) sum_rates = zeros(size(rates, 2)) @@ -337,7 +337,7 @@ function HopRatesGraphDsLij(species_hop_constants::Vector{F}, end function sample_target_site(hop_rates::HopRatesGraphDsLij, site, species, rng, - spatial_system) + spatial_system) @inbounds cum_hop_consts = hop_rates.hop_const_cumulative_sums[site] @inbounds n = searchsortedfirst(cum_hop_consts, rand(rng) * cum_hop_consts[end]) return nth_nbr(spatial_system, site, n) @@ -366,13 +366,13 @@ end function Base.show(io::IO, ::MIME"text/plain", hop_rates::HopRatesGridDsLij) num_specs, num_sites = length(hop_rates.species_hop_constants), - size(hop_rates.hop_const_cumulative_sums, 2) + size(hop_rates.hop_const_cumulative_sums, 2) println(io, - "HopRates with $num_specs species and $num_sites sites, optimized for CartesianGrid. \nHopping constants of form D_s * L_{i,j} where s is species, i is source and j is destination.") + "HopRates with $num_specs species and $num_sites sites, optimized for CartesianGrid. \nHopping constants of form D_s * L_{i,j} where s is species, i is source and j is destination.") end function HopRatesGridDsLij(species_hop_constants::Vector{F}, site_hop_constants::Matrix{F}; - do_cumsum = true) where {F <: Number} + do_cumsum = true) where {F <: Number} do_cumsum && (site_hop_constants = mapslices(cumsum, site_hop_constants, dims = 1)) rates = zeros(F, length(species_hop_constants), size(site_hop_constants, 2)) sum_rates = zeros(size(rates, 2)) @@ -380,9 +380,9 @@ function HopRatesGridDsLij(species_hop_constants::Vector{F}, site_hop_constants: end function HopRatesGridDsLij(species_hop_constants::Vector{F}, - site_hop_constants::Vector{Vector{F}}, grid) where {F <: Number} + site_hop_constants::Vector{Vector{F}}, grid) where {F <: Number} new_hopping_constants = Matrix{F}(undef, 2 * dimension(grid), - length(site_hop_constants)) + length(site_hop_constants)) for site in 1:length(site_hop_constants) nb_constants = @view new_hopping_constants[:, site] pad_hop_vec!(nb_constants, grid, site, site_hop_constants[site]) @@ -420,12 +420,12 @@ end function Base.show(io::IO, ::MIME"text/plain", hop_rates::HopRatesGraphDsiLij) num_specs, num_sites = size(hop_rates.species_hop_constants) println(io, - "HopRates with $num_specs species and $num_sites sites. \nHopping constants of form D_{s,i} * L_{i,j} where s is species, i is source and j is destination.") + "HopRates with $num_specs species and $num_sites sites. \nHopping constants of form D_{s,i} * L_{i,j} where s is species, i is source and j is destination.") end function HopRatesGraphDsiLij(species_hop_constants::Matrix{F}, - site_hop_constants::Vector{Vector{F}}; - do_cumsum = true) where {F <: Number} + site_hop_constants::Vector{Vector{F}}; + do_cumsum = true) where {F <: Number} @assert size(species_hop_constants, 2) == length(site_hop_constants) do_cumsum && (site_hop_constants = map(cumsum, site_hop_constants)) rates = zeros(F, length(species_hop_constants), length(site_hop_constants)) @@ -434,7 +434,7 @@ function HopRatesGraphDsiLij(species_hop_constants::Matrix{F}, end function sample_target_site(hop_rates::HopRatesGraphDsiLij, site, species, rng, - spatial_system) + spatial_system) @inbounds cum_hop_consts = hop_rates.hop_const_cumulative_sums[site] @inbounds n = searchsortedfirst(cum_hop_consts, rand(rng) * cum_hop_consts[end]) return nth_nbr(spatial_system, site, n) @@ -463,13 +463,14 @@ end function Base.show(io::IO, ::MIME"text/plain", hop_rates::HopRatesGridDsiLij) num_specs, num_sites = length(hop_rates.species_hop_constants), - size(hop_rates.hop_const_cumulative_sums, 2) + size(hop_rates.hop_const_cumulative_sums, 2) println(io, - "HopRates with $num_specs species and $num_sites sites, optimized for CartesianGrid. \nHopping constants of form D_{s,i} * L_{i,j} where s is species, i is source and j is destination.") + "HopRates with $num_specs species and $num_sites sites, optimized for CartesianGrid. \nHopping constants of form D_{s,i} * L_{i,j} where s is species, i is source and j is destination.") end -function HopRatesGridDsiLij(species_hop_constants::Matrix{F}, site_hop_constants::Matrix{F}; - do_cumsum = true) where {F <: Number} +function HopRatesGridDsiLij( + species_hop_constants::Matrix{F}, site_hop_constants::Matrix{F}; + do_cumsum = true) where {F <: Number} @assert size(species_hop_constants, 2) == size(site_hop_constants, 2) do_cumsum && (site_hop_constants = mapslices(cumsum, site_hop_constants, dims = 1)) rates = zeros(F, size(species_hop_constants)) @@ -478,9 +479,9 @@ function HopRatesGridDsiLij(species_hop_constants::Matrix{F}, site_hop_constants end function HopRatesGridDsiLij(species_hop_constants::Matrix{F}, - site_hop_constants::Vector{Vector{F}}, grid) where {F <: Number} + site_hop_constants::Vector{Vector{F}}, grid) where {F <: Number} new_hopping_constants = Matrix{F}(undef, 2 * dimension(grid), - length(site_hop_constants)) + length(site_hop_constants)) for site in 1:length(site_hop_constants) nb_constants = @view new_hopping_constants[:, site] pad_hop_vec!(nb_constants, grid, site, site_hop_constants[site]) diff --git a/src/spatial/nsm.jl b/src/spatial/nsm.jl index 422c2bf1..85e5cc42 100644 --- a/src/spatial/nsm.jl +++ b/src/spatial/nsm.jl @@ -4,7 +4,7 @@ #NOTE state vector u is a matrix. u[i,j] is species i, site j #NOTE hopping_constants is a matrix. hopping_constants[i,j] is species i, site j mutable struct NSMJumpAggregation{T, S, F1, F2, RNG, J, RX, HOP, DEPGR, VJMAP, JVMAP, - PQ, SS} <: + PQ, SS} <: AbstractSSAJumpAggregator{T, S, F1, F2, RNG} next_jump::SpatialJump{J} #some structure to identify the next event: reaction or hop prev_jump::SpatialJump{J} #some structure to identify the previous event: reaction or hop @@ -24,11 +24,12 @@ mutable struct NSMJumpAggregation{T, S, F1, F2, RNG, J, RX, HOP, DEPGR, VJMAP, J numspecies::Int #number of species end -function NSMJumpAggregation(nj::SpatialJump{J}, njt::T, et::T, rx_rates::RX, hop_rates::HOP, - sps::Tuple{Bool, Bool}, - rng::RNG, spatial_system::SS; num_specs, - vartojumps_map = nothing, jumptovars_map = nothing, - dep_graph = nothing, kwargs...) where {J, T, RX, HOP, RNG, SS} +function NSMJumpAggregation( + nj::SpatialJump{J}, njt::T, et::T, rx_rates::RX, hop_rates::HOP, + sps::Tuple{Bool, Bool}, + rng::RNG, spatial_system::SS; num_specs, + vartojumps_map = nothing, jumptovars_map = nothing, + dep_graph = nothing, kwargs...) where {J, T, RX, HOP, RNG, SS} # a dependency graph is needed if dep_graph === nothing @@ -55,30 +56,30 @@ function NSMJumpAggregation(nj::SpatialJump{J}, njt::T, et::T, rx_rates::RX, hop pq = MutableBinaryMinHeap{T}() NSMJumpAggregation{T, Nothing, Nothing, Nothing, RNG, J, RX, HOP, typeof(dg), - typeof(vtoj_map), typeof(jtov_map), typeof(pq), SS}(nj, nj, njt, et, - rx_rates, - hop_rates, - nothing, - nothing, - sps, rng, dg, - vtoj_map, - jtov_map, - pq, - spatial_system, - num_specs) + typeof(vtoj_map), typeof(jtov_map), typeof(pq), SS}(nj, nj, njt, et, + rx_rates, + hop_rates, + nothing, + nothing, + sps, rng, dg, + vtoj_map, + jtov_map, + pq, + spatial_system, + num_specs) end ############################# Required Functions ############################## # creating the JumpAggregation structure (function wrapper-based constant jumps) function aggregate(aggregator::NSM, starting_state, p, t, end_time, constant_jumps, - ma_jumps, save_positions, rng; hopping_constants, spatial_system, - kwargs...) + ma_jumps, save_positions, rng; hopping_constants, spatial_system, + kwargs...) num_species = size(starting_state, 1) majumps = ma_jumps if majumps === nothing majumps = MassActionJump(Vector{typeof(end_time)}(), - Vector{Vector{Pair{Int, Int}}}(), - Vector{Vector{Pair{Int, Int}}}()) + Vector{Vector{Pair{Int, Int}}}(), + Vector{Vector{Pair{Int, Int}}}()) end next_jump = SpatialJump{Int}(typemax(Int), typemax(Int), typemax(Int)) #a placeholder @@ -87,8 +88,8 @@ function aggregate(aggregator::NSM, starting_state, p, t, end_time, constant_jum hop_rates = HopRates(hopping_constants, spatial_system) NSMJumpAggregation(next_jump, next_jump_time, end_time, rx_rates, hop_rates, - save_positions, rng, spatial_system; num_specs = num_species, - kwargs...) + save_positions, rng, spatial_system; num_specs = num_species, + kwargs...) end # set up a new simulation and calculate the first jump / jump time diff --git a/src/spatial/reaction_rates.jl b/src/spatial/reaction_rates.jl index 8084166c..737cc5c9 100644 --- a/src/spatial/reaction_rates.jl +++ b/src/spatial/reaction_rates.jl @@ -53,7 +53,7 @@ end update rates of all reactions in rxs at site """ function update_rx_rates!(rx_rates::RxRates, rxs, u::AbstractMatrix, integrator, - site) + site) ma_jumps = rx_rates.ma_jumps @inbounds for rx in rxs rate = eval_massaction_rate(u, rx, ma_jumps, site) @@ -62,7 +62,7 @@ function update_rx_rates!(rx_rates::RxRates, rxs, u::AbstractMatrix, integrator, end function update_rx_rates!(rx_rates::RxRates, rxs, integrator, - site) + site) u = integrator.u update_rx_rates!(rx_rates, rxs, u, integrator, site) end @@ -74,7 +74,7 @@ sample a reaction at site, return reaction index """ function sample_rx_at_site(rx_rates::RxRates, site, rng) linear_search((@view rx_rates.rates[:, site]), - rand(rng) * total_site_rx_rate(rx_rates, site)) + rand(rng) * total_site_rx_rate(rx_rates, site)) end # helper functions diff --git a/src/spatial/spatial_massaction_jump.jl b/src/spatial/spatial_massaction_jump.jl index 4379804b..f1a4fea5 100644 --- a/src/spatial/spatial_massaction_jump.jl +++ b/src/spatial/spatial_massaction_jump.jl @@ -13,12 +13,12 @@ struct SpatialMassActionJump{A <: AVecOrNothing, B <: AMatOrNothing, S, U, V} <: uniform rates go first in ordering """ function SpatialMassActionJump{A, B, S, U, V}(uniform_rates::A, spatial_rates::B, - reactant_stoch::S, net_stoch::U, - param_mapper::V, scale_rates::Bool, - useiszero::Bool, - nocopy::Bool) where {A <: AVecOrNothing, - B <: AMatOrNothing, - S, U, V} + reactant_stoch::S, net_stoch::U, + param_mapper::V, scale_rates::Bool, + useiszero::Bool, + nocopy::Bool) where {A <: AVecOrNothing, + B <: AMatOrNothing, + S, U, V} uniform_rates = (nocopy || isnothing(uniform_rates)) ? uniform_rates : copy(uniform_rates) spatial_rates = (nocopy || isnothing(spatial_rates)) ? spatial_rates : @@ -44,93 +44,93 @@ end ################ Constructors ################## function SpatialMassActionJump(urates::A, srates::B, rs::S, ns::U, pmapper::V; - scale_rates = true, useiszero = true, - nocopy = false) where {A <: AVecOrNothing, - B <: AMatOrNothing, S, U, V} + scale_rates = true, useiszero = true, + nocopy = false) where {A <: AVecOrNothing, + B <: AMatOrNothing, S, U, V} SpatialMassActionJump{A, B, S, U, V}(urates, srates, rs, ns, pmapper, scale_rates, - useiszero, nocopy) + useiszero, nocopy) end function SpatialMassActionJump(urates::A, srates::B, rs, ns; scale_rates = true, - useiszero = true, - nocopy = false) where {A <: AVecOrNothing, - B <: AMatOrNothing} + useiszero = true, + nocopy = false) where {A <: AVecOrNothing, + B <: AMatOrNothing} SpatialMassActionJump(urates, srates, rs, ns, nothing; scale_rates = scale_rates, - useiszero = useiszero, nocopy = nocopy) + useiszero = useiszero, nocopy = nocopy) end function SpatialMassActionJump(srates::B, rs, ns, pmapper; scale_rates = true, - useiszero = true, - nocopy = false) where {B <: AMatOrNothing} + useiszero = true, + nocopy = false) where {B <: AMatOrNothing} SpatialMassActionJump(nothing, srates, rs, ns, pmapper; scale_rates = scale_rates, - useiszero = useiszero, nocopy = nocopy) + useiszero = useiszero, nocopy = nocopy) end function SpatialMassActionJump(srates::B, rs, ns; scale_rates = true, useiszero = true, - nocopy = false) where {B <: AMatOrNothing} + nocopy = false) where {B <: AMatOrNothing} SpatialMassActionJump(nothing, srates, rs, ns, nothing; scale_rates = scale_rates, - useiszero = useiszero, nocopy = nocopy) + useiszero = useiszero, nocopy = nocopy) end function SpatialMassActionJump(urates::A, rs, ns, pmapper; scale_rates = true, - useiszero = true, - nocopy = false) where {A <: AVecOrNothing} + useiszero = true, + nocopy = false) where {A <: AVecOrNothing} SpatialMassActionJump(urates, nothing, rs, ns, pmapper; scale_rates = scale_rates, - useiszero = useiszero, nocopy = nocopy) + useiszero = useiszero, nocopy = nocopy) end function SpatialMassActionJump(urates::A, rs, ns; scale_rates = true, useiszero = true, - nocopy = false) where {A <: AVecOrNothing} + nocopy = false) where {A <: AVecOrNothing} SpatialMassActionJump(urates, nothing, rs, ns, nothing; scale_rates = scale_rates, - useiszero = useiszero, nocopy = nocopy) + useiszero = useiszero, nocopy = nocopy) end function SpatialMassActionJump(ma_jumps::MassActionJump{T, S, U, V}; scale_rates = true, - useiszero = true, nocopy = false) where {T, S, U, V} + useiszero = true, nocopy = false) where {T, S, U, V} SpatialMassActionJump(ma_jumps.scaled_rates, ma_jumps.reactant_stoch, - ma_jumps.net_stoch, ma_jumps.param_mapper; - scale_rates = scale_rates, useiszero = useiszero, nocopy = nocopy) + ma_jumps.net_stoch, ma_jumps.param_mapper; + scale_rates = scale_rates, useiszero = useiszero, nocopy = nocopy) end ############################################## function get_num_majumps(smaj::SpatialMassActionJump{ - Nothing, Nothing, S, U, V}) where - {S, U, V} + Nothing, Nothing, S, U, V}) where + {S, U, V} 0 end function get_num_majumps(smaj::SpatialMassActionJump{ - Nothing, B, S, U, V}) where - {B, S, U, V} + Nothing, B, S, U, V}) where + {B, S, U, V} size(smaj.spatial_rates, 1) end function get_num_majumps(smaj::SpatialMassActionJump{ - A, Nothing, S, U, V}) where - {A, S, U, V} + A, Nothing, S, U, V}) where + {A, S, U, V} length(smaj.uniform_rates) end function get_num_majumps(smaj::SpatialMassActionJump{ - A, B, S, U, V}) where - {A <: AbstractVector, B <: AbstractMatrix, S, U, V} + A, B, S, U, V}) where + {A <: AbstractVector, B <: AbstractMatrix, S, U, V} length(smaj.uniform_rates) + size(smaj.spatial_rates, 1) end using_params(smaj::SpatialMassActionJump) = false function rate_at_site(rx, site, - smaj::SpatialMassActionJump{Nothing, B, S, U, V}) where {B, S, U, V} + smaj::SpatialMassActionJump{Nothing, B, S, U, V}) where {B, S, U, V} smaj.spatial_rates[rx, site] end function rate_at_site(rx, site, - smaj::SpatialMassActionJump{A, Nothing, S, U, V}) where {A, S, U, V} + smaj::SpatialMassActionJump{A, Nothing, S, U, V}) where {A, S, U, V} smaj.uniform_rates[rx] end function rate_at_site(rx, site, - smaj::SpatialMassActionJump{A, B, S, U, V}) where - {A <: AbstractVector, B <: AbstractMatrix, S, U, V} + smaj::SpatialMassActionJump{A, B, S, U, V}) where + {A <: AbstractVector, B <: AbstractMatrix, S, U, V} num_unif_rxs = length(smaj.uniform_rates) rx <= num_unif_rxs ? smaj.uniform_rates[rx] : smaj.spatial_rates[rx - num_unif_rxs, site] end function evalrxrate(speciesmat::AbstractMatrix{T}, rxidx::S, majump::SpatialMassActionJump, - site::Int) where {T, S} + site::Int) where {T, S} val = one(T) @inbounds for specstoch in majump.reactant_stoch[rxidx] specpop = speciesmat[specstoch[1], site] diff --git a/src/spatial/topology.jl b/src/spatial/topology.jl index 3bb0b901..9ad6381e 100644 --- a/src/spatial/topology.jl +++ b/src/spatial/topology.jl @@ -15,7 +15,7 @@ const offsets_2D = [ CartesianIndex(0, -1), CartesianIndex(-1, 0), CartesianIndex(1, 0), - CartesianIndex(0, 1), + CartesianIndex(0, 1) ] const offsets_3D = [ CartesianIndex(0, 0, -1), @@ -23,7 +23,7 @@ const offsets_3D = [ CartesianIndex(-1, 0, 0), CartesianIndex(1, 0, 0), CartesianIndex(0, 1, 0), - CartesianIndex(0, 0, 1), + CartesianIndex(0, 0, 1) ] """ @@ -86,7 +86,7 @@ end given a vector of hopping constants of length num_neighbors(grid, site), form the vector of size 2*(dimension of grid) with zeros at indices where the neighbor is out of bounds. Store it in to_pad """ function pad_hop_vec!(to_pad::AbstractVector{F}, grid, site, - hop_vec::Vector{F}) where {F <: Number} + hop_vec::Vector{F}) where {F <: Number} CI = grid.CI I = CI[site] nbr_counter = 1 @@ -145,6 +145,6 @@ function rand_nbr(rng, grid::CartesianGridRej, site::Int) end function Base.show(io::IO, ::MIME"text/plain", - grid::CartesianGridRej) + grid::CartesianGridRej) println(io, "A Cartesian grid with dimensions $(grid.dims)") end diff --git a/src/spatial/utils.jl b/src/spatial/utils.jl index d153dd54..370e8601 100644 --- a/src/spatial/utils.jl +++ b/src/spatial/utils.jl @@ -18,7 +18,7 @@ end function Base.show(io::IO, ::MIME"text/plain", jump::SpatialJump) println(io, - "SpatialJump with source $(jump.src), destination $(jump.dst) and index $(jump.jidx).") + "SpatialJump with source $(jump.src), destination $(jump.dst) and index $(jump.jidx).") end ######################## helper routines for all spatial SSAs ######################## @@ -35,7 +35,7 @@ function sample_jump_direct(p, site) else species_to_diffuse, target_site = sample_hop_at_site(p.hop_rates, site, p.rng, - p.spatial_system) + p.spatial_system) return SpatialJump(site, species_to_diffuse, target_site) end end @@ -71,7 +71,7 @@ function update_state!(p, integrator) else rx_index = reaction_id_from_jump(p, jump) @inbounds executerx!((@view integrator.u[:, jump.src]), rx_index, - p.rx_rates.ma_jumps) + p.rx_rates.ma_jumps) end # save jump that was just executed p.prev_jump = jump diff --git a/test/allocations.jl b/test/allocations.jl index 9bf41a65..64efb758 100644 --- a/test/allocations.jl +++ b/test/allocations.jl @@ -56,7 +56,7 @@ let end function makeprob(; T = 100.0, alg = Direct(), save_positions = (false, false), - graphkwargs = (;), rng) + graphkwargs = (;), rng) r1(u, p, t) = rate(p[1], u[1], u[2], p[2]) * u[1] r2(u, p, t) = rate(p[1], u[2], u[1], p[2]) * u[2] r3(u, p, t) = p[3] * u[1] @@ -82,11 +82,11 @@ let dprob = DiscreteProblem(u0, tspan, p) jprob = JumpProblem(dprob, alg, - ConstantRateJump(r1, aff1!), ConstantRateJump(r2, aff2!), - ConstantRateJump(r3, aff3!), - ConstantRateJump(r4, aff4!), ConstantRateJump(r5, aff5!), - ConstantRateJump(r6, aff6!); - save_positions, rng, graphkwargs...) + ConstantRateJump(r1, aff1!), ConstantRateJump(r2, aff2!), + ConstantRateJump(r3, aff3!), + ConstantRateJump(r4, aff4!), ConstantRateJump(r5, aff5!), + ConstantRateJump(r6, aff6!); + save_positions, rng, graphkwargs...) return jprob end diff --git a/test/bimolerx_test.jl b/test/bimolerx_test.jl index c6021dac..befeefab 100644 --- a/test/bimolerx_test.jl +++ b/test/bimolerx_test.jl @@ -40,14 +40,14 @@ reactstoch = [ [2 => 1], [1 => 1, 2 => 1], [3 => 1], - [3 => 3], + [3 => 3] ] netstoch = [ [1 => -2, 2 => 1], [1 => 2, 2 => -1], [1 => -1, 2 => -1, 3 => 1], [1 => 1, 2 => 1, 3 => -1], - [1 => 3, 3 => -3], + [1 => 3, 3 => -3] ] rates = [1.0, 2.0, 0.5, 0.75, 0.25] spec_to_dep_jumps = [[1, 3], [2, 3], [4, 5]] @@ -71,8 +71,8 @@ prob = DiscreteProblem(u0, (0.0, tf), rates) if doplot for alg in SSAalgs local jump_prob = JumpProblem(prob, alg, majumps, - vartojumps_map = spec_to_dep_jumps, - jumptovars_map = jump_to_dep_specs, rng = rng) + vartojumps_map = spec_to_dep_jumps, + jumptovars_map = jump_to_dep_specs, rng = rng) local sol = solve(jump_prob, SSAStepper()) local plothand = plot(sol, seriestype = :steppost, reuse = false) display(plothand) @@ -83,12 +83,12 @@ end if dotestmean for (i, alg) in enumerate(SSAalgs) local jump_prob = JumpProblem(prob, alg, majumps, save_positions = (false, false), - vartojumps_map = spec_to_dep_jumps, - jumptovars_map = jump_to_dep_specs, rng = rng) + vartojumps_map = spec_to_dep_jumps, + jumptovars_map = jump_to_dep_specs, rng = rng) means = runSSAs(jump_prob) relerr = abs(means - expected_avg) / expected_avg doprintmeans && println("Mean from method: ", typeof(alg), " is = ", means, - ", rel err = ", relerr) + ", rel err = ", relerr) @test abs(means - expected_avg) < reltol * expected_avg # test not specifying SSAStepper @@ -106,13 +106,13 @@ if dotestmean end jset = JumpSet((), (), nothing, majump_vec) jump_prob = JumpProblem(prob, Direct(), jset, save_positions = (false, false), - vartojumps_map = spec_to_dep_jumps, - jumptovars_map = jump_to_dep_specs, rng = rng) + vartojumps_map = spec_to_dep_jumps, + jumptovars_map = jump_to_dep_specs, rng = rng) meanval = runSSAs(jump_prob) relerr = abs(meanval - expected_avg) / expected_avg if doprintmeans println("Using individual MassActionJumps; Mean from method: ", typeof(Direct()), - " is = ", meanval, ", rel err = ", relerr) + " is = ", meanval, ", rel err = ", relerr) end @test abs(meanval - expected_avg) < reltol * expected_avg end diff --git a/test/bracketing.jl b/test/bracketing.jl index 68021457..66ba1cbe 100644 --- a/test/bracketing.jl +++ b/test/bracketing.jl @@ -35,7 +35,7 @@ majump_rates = [0.1] # death at rate 0.1 reactstoch = [[1 => 1]] netstoch = [[1 => -1]] majump = MassActionJump(majump_rates, reactstoch, - netstoch) + netstoch) reaction_index = 1 @test JP.get_majump_brackets(ulow, uhigh, reaction_index, majump)[1] == majump_rates[1] * ulow[1] # low diff --git a/test/degenerate_rx_cases.jl b/test/degenerate_rx_cases.jl index 97eda232..de6938c0 100644 --- a/test/degenerate_rx_cases.jl +++ b/test/degenerate_rx_cases.jl @@ -93,12 +93,12 @@ end dep_graph = [ [1, 2], - [1, 2], + [1, 2] ] spec_to_dep_jumps = [[2]] jump_to_dep_specs = [[1], [1]] namedpars = (dep_graph = dep_graph, vartojumps_map = spec_to_dep_jumps, - jumptovars_map = jump_to_dep_specs) + jumptovars_map = jump_to_dep_specs) for method in methods local jump_prob = JumpProblem(prob, method, jump, jump2; rng = rng, namedpars...) @@ -110,7 +110,7 @@ for method in methods if doprint println("Mix of constant and mass action jumps, method = ", typeof(method), - ", sol[end] = ", sol[end, end]) + ", sol[end] = ", sol[end, end]) end @test sol[end, end] > 200 end @@ -131,11 +131,11 @@ maj2 = MassActionJump(rs2, ns2; param_idxs = 2) js = JumpSet(maj1, maj2) maj = MassActionJump([rs1, rs2], [ns1, ns2]; param_idxs = [1, 2]) @test all(getfield(maj, fn) == getfield(js.massaction_jump, fn) - for fn in [:scaled_rates, :reactant_stoch, :net_stoch]) +for fn in [:scaled_rates, :reactant_stoch, :net_stoch]) @test all(maj.param_mapper.param_idxs .== js.massaction_jump.param_mapper.param_idxs) maj1 = MassActionJump([rs1], [ns1]; param_idxs = [1]) maj2 = MassActionJump([rs2], [ns2]; param_idxs = [2]) js = JumpSet(maj1, maj2) @test all(getfield(maj, fn) == getfield(js.massaction_jump, fn) - for fn in [:scaled_rates, :reactant_stoch, :net_stoch]) +for fn in [:scaled_rates, :reactant_stoch, :net_stoch]) @test all(maj.param_mapper.param_idxs .== js.massaction_jump.param_mapper.param_idxs) diff --git a/test/extended_jump_array.jl b/test/extended_jump_array.jl index 553ee514..466ab690 100644 --- a/test/extended_jump_array.jl +++ b/test/extended_jump_array.jl @@ -5,7 +5,7 @@ rng = StableRNG(123) # Check that the new broadcast norm gives the same result as the old one rand_array = ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}}(rand(rng, 5), - rand(rng, 2)) + rand(rng, 2)) old_norm = Base.FastMath.sqrt_fast(DiffEqBase.UNITLESS_ABS2(rand_array) / max(DiffEqBase.recursive_length(rand_array), 1)) new_norm = DiffEqBase.ODE_DEFAULT_NORM(rand_array, 0.0) @@ -13,8 +13,8 @@ new_norm = DiffEqBase.ODE_DEFAULT_NORM(rand_array, 0.0) # Check for an ExtendedJumpArray where the types differ (Float64/Int64) rand_array = ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Int64}}(rand(rng, 5), - rand(rng, 1:1000, - 2)) + rand(rng, 1:1000, + 2)) old_norm = Base.FastMath.sqrt_fast(DiffEqBase.UNITLESS_ABS2(rand_array) / max(DiffEqBase.recursive_length(rand_array), 1)) new_norm = DiffEqBase.ODE_DEFAULT_NORM(rand_array, 0.0) diff --git a/test/extinction_test.jl b/test/extinction_test.jl index 857da667..880254ba 100644 --- a/test/extinction_test.jl +++ b/test/extinction_test.jl @@ -22,7 +22,7 @@ algs = (JumpProcesses.JUMP_AGGREGATORS..., JumpProcesses.NullAggregator()) for n in 1:Nsims for ssa in algs local jprob = JumpProblem(dprob, ssa, majump, save_positions = (false, false), - rng = rng) + rng = rng) local sol = solve(jprob, SSAStepper()) @test sol[1, end] == 0 @test sol.t[end] < Inf @@ -34,7 +34,7 @@ dprob = DiscreteProblem(u0, (0.0, 100.0), rates) for ssa in algs local jprob = JumpProblem(dprob, ssa, majump, save_positions = (false, false), - rng = rng) + rng = rng) local sol = solve(jprob, SSAStepper(), saveat = 100.0) @test sol[1, end] == 0 @test sol.t[end] < Inf @@ -71,7 +71,7 @@ function extinction_affect!2(integrator) nothing end cb = DiscreteCallback(extinction_condition2, extinction_affect!2, - save_positions = (false, false)) + save_positions = (false, false)) dprob = DiscreteProblem(u0, (0.0, 1000.0), rates) jprob = JumpProblem(dprob, majump; save_positions = (false, false), rng) sol = solve(jprob; callback = cb, save_end = false) diff --git a/test/functionwrappers.jl b/test/functionwrappers.jl index d2d0f2b4..2f009ead 100644 --- a/test/functionwrappers.jl +++ b/test/functionwrappers.jl @@ -9,7 +9,7 @@ let nothing end jump = VariableRateJump(rate, affect!; urate = (u, p, t) -> 10.0, - rateinterval = (u, p, t) -> 0.1) + rateinterval = (u, p, t) -> 0.1) prob = DiscreteProblem([0.0], (0.0, 2.0), [1.0]) jprob = JumpProblem(prob, Coevolve(), jump; dep_graph = [[1]], rng) diff --git a/test/geneexpr_test.jl b/test/geneexpr_test.jl index 3cc65915..07a49bba 100644 --- a/test/geneexpr_test.jl +++ b/test/geneexpr_test.jl @@ -61,7 +61,7 @@ reactstoch = [ [2 => 1], [3 => 1], [1 => 1, 3 => 1], - [4 => 1], + [4 => 1] ] netstoch = [ [2 => 1], @@ -69,7 +69,7 @@ netstoch = [ [2 => -1], [3 => -1], [1 => -1, 3 => -1, 4 => 1], - [1 => 1, 3 => 1, 4 => -1], + [1 => 1, 3 => 1, 4 => -1] ] spec_to_dep_jumps = [[1, 5], [2, 3], [4, 5], [6]] jump_to_dep_specs = [[2], [3], [2], [3], [1, 3, 4], [1, 3, 4]] @@ -88,8 +88,8 @@ if doplot plothand = plot(reuse = false) for alg in SSAalgs local jump_prob = JumpProblem(prob, alg, majumps, - vartojumps_map = spec_to_dep_jumps, - jumptovars_map = jump_to_dep_specs, rng = rng) + vartojumps_map = spec_to_dep_jumps, + jumptovars_map = jump_to_dep_specs, rng = rng) local sol = solve(jump_prob, SSAStepper()) plot!(plothand, sol.t, sol[3, :], seriestype = :steppost) end @@ -100,12 +100,12 @@ end if dotestmean for (i, alg) in enumerate(SSAalgs) local jump_prob = JumpProblem(prob, alg, majumps, save_positions = (false, false), - vartojumps_map = spec_to_dep_jumps, - jumptovars_map = jump_to_dep_specs, rng = rng) + vartojumps_map = spec_to_dep_jumps, + jumptovars_map = jump_to_dep_specs, rng = rng) means = runSSAs(jump_prob) relerr = abs(means - expected_avg) / expected_avg doprintmeans && println("Mean from method: ", typeof(alg), " is = ", means, - ", rel err = ", relerr) + ", rel err = ", relerr) @test abs(means - expected_avg) < reltol * expected_avg means = runSSAs(jump_prob; use_stepper = false) @@ -113,20 +113,20 @@ if dotestmean @test abs(means - expected_avg) < reltol * expected_avg jump_probf = JumpProblem(probf, alg, majumps, save_positions = (false, false), - vartojumps_map = spec_to_dep_jumps, - jumptovars_map = jump_to_dep_specs, rng = rng) + vartojumps_map = spec_to_dep_jumps, + jumptovars_map = jump_to_dep_specs, rng = rng) means = runSSAs(jump_probf) relerr = abs(means - expected_avg) / expected_avg doprintmeans && println("Mean from method: ", typeof(alg), " is = ", means, - ", rel err = ", relerr) + ", rel err = ", relerr) @test abs(means - expected_avg) < reltol * expected_avg end end # no-aggregator tests jump_prob = JumpProblem(prob, majumps; save_positions = (false, false), - vartojumps_map = spec_to_dep_jumps, - jumptovars_map = jump_to_dep_specs, rng) + vartojumps_map = spec_to_dep_jumps, + jumptovars_map = jump_to_dep_specs, rng) @test abs(runSSAs(jump_prob) - expected_avg) < reltol * expected_avg @test abs(runSSAs(jump_prob; use_stepper = false) - expected_avg) < reltol * expected_avg @@ -166,15 +166,15 @@ let nothing end crjs = JumpSet(ConstantRateJump(r1, a1!), ConstantRateJump(r2, a2!), - ConstantRateJump(r3, a3!), ConstantRateJump(r4, a4!), - ConstantRateJump(r5, a5!), - ConstantRateJump(r6, a6!)) + ConstantRateJump(r3, a3!), ConstantRateJump(r4, a4!), + ConstantRateJump(r5, a5!), + ConstantRateJump(r6, a6!)) vrjs = JumpSet(VariableRateJump(r1, a1!; save_positions = (false, false)), - VariableRateJump(r2, a2!, save_positions = (false, false)), - VariableRateJump(r3, a3!, save_positions = (false, false)), - VariableRateJump(r4, a4!, save_positions = (false, false)), - VariableRateJump(r5, a5!, save_positions = (false, false)), - VariableRateJump(r6, a6!, save_positions = (false, false))) + VariableRateJump(r2, a2!, save_positions = (false, false)), + VariableRateJump(r3, a3!, save_positions = (false, false)), + VariableRateJump(r4, a4!, save_positions = (false, false)), + VariableRateJump(r5, a5!, save_positions = (false, false)), + VariableRateJump(r6, a6!, save_positions = (false, false))) prob = DiscreteProblem(u0, (0.0, tf), rates) crjprob = JumpProblem(prob, crjs; save_positions = (false, false), rng) diff --git a/test/hawkes_test.jl b/test/hawkes_test.jl index 977f1665..7e4623dd 100644 --- a/test/hawkes_test.jl +++ b/test/hawkes_test.jl @@ -61,9 +61,9 @@ function hawkes_jump(u, g, h; uselrate = true) end function hawkes_problem(p, agg::Coevolve; u = [0.0], - tspan = (0.0, 50.0), - save_positions = (false, true), - g = [[1]], h = [[]], uselrate = true) + tspan = (0.0, 50.0), + save_positions = (false, true), + g = [[1]], h = [[]], uselrate = true) dprob = DiscreteProblem(u, tspan, p) jumps = hawkes_jump(u, g, h; uselrate) jprob = JumpProblem(dprob, agg, jumps...; dep_graph = g, save_positions, rng) @@ -76,8 +76,8 @@ function f!(du, u, p, t) end function hawkes_problem(p, agg; u = [0.0], tspan = (0.0, 50.0), - save_positions = (false, true), - g = [[1]], h = [[]], kwargs...) + save_positions = (false, true), + g = [[1]], h = [[]], kwargs...) oprob = ODEProblem(f!, u, tspan, p) jumps = hawkes_jump(u, g, h) jprob = JumpProblem(oprob, agg, jumps...; save_positions, rng) @@ -154,7 +154,7 @@ let alg = Coevolve() oprob = ODEProblem(f!, u0, tspan, p) jumps = hawkes_jump(u0, g, h) jprob = JumpProblem(oprob, alg, jumps...; dep_graph = g, rng, - use_vrj_bounds = false) + use_vrj_bounds = false) @test length(jprob.variable_jumps) == 1 sols = Vector{ODESolution}(undef, Nsims) for n in 1:Nsims diff --git a/test/jprob_symbol_indexing.jl b/test/jprob_symbol_indexing.jl index 2915f939..4cad845c 100644 --- a/test/jprob_symbol_indexing.jl +++ b/test/jprob_symbol_indexing.jl @@ -8,7 +8,7 @@ crj1 = ConstantRateJump(rate1, affect1!) crj2 = ConstantRateJump(rate2, affect2!) maj = MassActionJump([[1 => 1], [1 => 1]], [[1 => -1], [1 => -1]]; param_idxs = [1, 2]) g = DiscreteFunction((du, u, p, t) -> nothing; - sys = SymbolicIndexingInterface.SymbolCache([:a, :b], [:p1, :p2], :t)) + sys = SymbolicIndexingInterface.SymbolCache([:a, :b], [:p1, :p2], :t)) dprob = DiscreteProblem(g, [0, 10], (0.0, 10.0), [1.0, 2.0]) jprob = JumpProblem(dprob, Direct(), crj1, crj2, maj) diff --git a/test/linearreaction_test.jl b/test/linearreaction_test.jl index 15ae0722..657982a6 100644 --- a/test/linearreaction_test.jl +++ b/test/linearreaction_test.jl @@ -53,7 +53,7 @@ function A_to_B_tuple(N, method) jset = JumpSet((), jumps, nothing, nothing) prob = DiscreteProblem([A0, 0], (0.0, tf)) jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng, - namedpars...) + namedpars...) jump_prob end @@ -75,7 +75,7 @@ function A_to_B_vec(N, method) jset = JumpSet((), jumps, nothing, nothing) prob = DiscreteProblem([A0, 0], (0.0, tf)) jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng, - namedpars...) + namedpars...) jump_prob end @@ -93,7 +93,7 @@ function A_to_B_ma(N, method) jset = JumpSet((), (), nothing, majumps) prob = DiscreteProblem([A0, 0], (0.0, tf)) jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng, - namedpars...) + namedpars...) jump_prob end @@ -127,7 +127,7 @@ function A_to_B_hybrid(N, method) jset = JumpSet((), jumps, nothing, majumps) prob = DiscreteProblem([A0, 0], (0.0, tf)) jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng, - namedpars...) + namedpars...) jump_prob end @@ -161,7 +161,7 @@ function A_to_B_hybrid_nojset(N, method) jumps = (constjumps..., majumps) prob = DiscreteProblem([A0, 0], (0.0, tf)) jump_prob = JumpProblem(prob, method, jumps...; save_positions = (false, false), - rng, namedpars...) + rng, namedpars...) jump_prob end @@ -191,7 +191,7 @@ function A_to_B_hybrid_vecs(N, method) jset = JumpSet((), jumpvec, nothing, majumps) prob = DiscreteProblem([A0, 0], (0.0, tf)) jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng, - namedpars...) + namedpars...) jump_prob end @@ -221,7 +221,7 @@ function A_to_B_hybrid_vecs_scalars(N, method) jset = JumpSet((), jumpvec, nothing, majumps) prob = DiscreteProblem([A0, 0], (0.0, tf)) jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng, - namedpars...) + namedpars...) jump_prob end @@ -252,7 +252,7 @@ function A_to_B_hybrid_tups_scalars(N, method) jumps = ((maj for maj in majumpsv)..., (jump for jump in jumpvec)...) prob = DiscreteProblem([A0, 0], (0.0, tf)) jump_prob = JumpProblem(prob, method, jumps...; save_positions = (false, false), - rng, namedpars...) + rng, namedpars...) jump_prob end @@ -283,7 +283,7 @@ function A_to_B_hybrid_tups(N, method) jset = JumpSet((), jumps, nothing, majumps) prob = DiscreteProblem([A0, 0], (0.0, tf)) jump_prob = JumpProblem(prob, method, jset; save_positions = (false, false), rng, - namedpars...) + namedpars...) jump_prob end @@ -303,7 +303,7 @@ for method in SSAalgs meanval = runSSAs(jump_prob) if doprint println("Method: ", method, ", Jump input types: ", jump_prob_gen, - ", sample mean = ", meanval, ", actual mean = ", exactmeanval) + ", sample mean = ", meanval, ", actual mean = ", exactmeanval) end @test abs(meanval - exactmeanval) < 1.0 end @@ -320,7 +320,7 @@ for method in SSAalgs meanval = runSSAs(jump_prob) if doprint println("Method: ", method, ", Jump input types: ", jump_prob_gen, - ", sample mean = ", meanval, ", actual mean = ", exactmeanval) + ", sample mean = ", meanval, ", actual mean = ", exactmeanval) end @test abs(meanval - exactmeanval) < 1.0 end diff --git a/test/monte_carlo_test.jl b/test/monte_carlo_test.jl index d3bf4971..0563ea07 100644 --- a/test/monte_carlo_test.jl +++ b/test/monte_carlo_test.jl @@ -11,12 +11,12 @@ jump = VariableRateJump(rate, affect!, save_positions = (false, true)) jump_prob = JumpProblem(prob, Direct(), jump; rng = rng) monte_prob = EnsembleProblem(jump_prob) sol = solve(monte_prob, SRIW1(), EnsembleSerial(), trajectories = 3, - save_everystep = false, dt = 0.001, adaptive = false) + save_everystep = false, dt = 0.001, adaptive = false) @test sol.u[1].t[2] != sol.u[2].t[2] != sol.u[3].t[2] jump = ConstantRateJump(rate, affect!) jump_prob = JumpProblem(prob, Direct(), jump, save_positions = (true, false), rng = rng) monte_prob = EnsembleProblem(jump_prob) sol = solve(monte_prob, SRIW1(), EnsembleSerial(), trajectories = 3, - save_everystep = false, dt = 0.001, adaptive = false) + save_everystep = false, dt = 0.001, adaptive = false) @test sol.u[1].t[2] != sol.u[2].t[2] != sol.u[3].t[2] diff --git a/test/remake_test.jl b/test/remake_test.jl index acb414a6..676e615d 100644 --- a/test/remake_test.jl +++ b/test/remake_test.jl @@ -22,7 +22,7 @@ tspan = (0.0, 2500.0) dprob = DiscreteProblem(u0, tspan, p) jprob = JumpProblem(dprob, Direct(), jump, jump2, save_positions = (false, false), - rng = rng) + rng = rng) sol = solve(jprob, SSAStepper()) @test sol[3, end] == 1000 diff --git a/test/reversible_binding.jl b/test/reversible_binding.jl index f019f6c4..f872d92a 100644 --- a/test/reversible_binding.jl +++ b/test/reversible_binding.jl @@ -8,11 +8,11 @@ Nsims = 1e4 # ABC model A + B <--> C reactstoch = [ [1 => 1, 2 => 1], - [3 => 1], + [3 => 1] ] netstoch = [ [1 => -1, 2 => -1, 3 => 1], - [1 => 1, 2 => 1, 3 => -1], + [1 => 1, 2 => 1, 3 => -1] ] rates = [0.1, 1.0] u0 = [500, 500, 0] @@ -49,7 +49,7 @@ algs = (JumpProcesses.JUMP_AGGREGATORS..., JumpProcesses.NullAggregator()) relative_tolerance = 0.01 for alg in algs local jprob = JumpProblem(prob, alg, majumps, save_positions = (false, false), - rng = rng) + rng = rng) local Amean = getmean(jprob, Nsims) @test abs(Amean - mastereq_mean) / mastereq_mean < relative_tolerance end diff --git a/test/save_positions.jl b/test/save_positions.jl index 450a1161..1e5ddc40 100644 --- a/test/save_positions.jl +++ b/test/save_positions.jl @@ -13,17 +13,17 @@ let # Coevolve will consider many candidates before the end of the simmulation. # None of these points should be saved. jump = VariableRateJump((u, p, t) -> 0, (integrator) -> integrator.u[1] += 1; - urate = (u, p, t) -> 1.0, rateinterval = (u, p, t) -> 5.0) + urate = (u, p, t) -> 1.0, rateinterval = (u, p, t) -> 5.0) jumpproblem = JumpProblem(dprob, alg, jump; dep_graph = [[1]], - save_positions = (false, true), rng) + save_positions = (false, true), rng) sol = solve(jumpproblem, SSAStepper()) @test sol.t == [0.0, 30.0] oprob = ODEProblem((du, u, p, t) -> 0, u0, tspan) jump = VariableRateJump((u, p, t) -> 0, (integrator) -> integrator.u[1] += 1; - urate = (u, p, t) -> 1.0, rateinterval = (u, p, t) -> 5.0) + urate = (u, p, t) -> 1.0, rateinterval = (u, p, t) -> 5.0) jumpproblem = JumpProblem(oprob, alg, jump; dep_graph = [[1]], - save_positions = (false, true), rng) + save_positions = (false, true), rng) sol = solve(jumpproblem, Tsit5(); save_everystep = false) @test sol.t == [0.0, 30.0] end diff --git a/test/saveat_regression.jl b/test/saveat_regression.jl index 4ee46cf7..03665d76 100644 --- a/test/saveat_regression.jl +++ b/test/saveat_regression.jl @@ -15,7 +15,7 @@ ts = collect(0:0.002:tspan[2]) NA = zeros(length(ts)) Nsims = 10_000 sol = JumpProcesses.solve(EnsembleProblem(jprob), SSAStepper(), saveat = ts, - trajectories = Nsims) + trajectories = Nsims) for i in 1:length(sol) NA .+= sol.u[i][1, :] diff --git a/test/spatial/ABC.jl b/test/spatial/ABC.jl index f01e1e61..8d7230a7 100644 --- a/test/spatial/ABC.jl +++ b/test/spatial/ABC.jl @@ -50,17 +50,17 @@ end # testing grids = [CartesianGridRej(dims), Graphs.grid(dims)] jump_problems = JumpProblem[JumpProblem(prob, NSM(), majumps, - hopping_constants = hopping_constants, - spatial_system = grid, - save_positions = (false, false), rng = rng) + hopping_constants = hopping_constants, + spatial_system = grid, + save_positions = (false, false), rng = rng) for grid in grids] push!(jump_problems, - JumpProblem(prob, DirectCRDirect(), majumps, hopping_constants = hopping_constants, - spatial_system = grids[1], save_positions = (false, false), rng = rng)) + JumpProblem(prob, DirectCRDirect(), majumps, hopping_constants = hopping_constants, + spatial_system = grids[1], save_positions = (false, false), rng = rng)) # setup flattenned jump prob push!(jump_problems, - JumpProblem(prob, NRM(), majumps, hopping_constants = hopping_constants, - spatial_system = grids[1], save_positions = (false, false), rng = rng)) + JumpProblem(prob, NRM(), majumps, hopping_constants = hopping_constants, + spatial_system = grids[1], save_positions = (false, false), rng = rng)) # test for spatial_jump_prob in jump_problems solution = solve(spatial_jump_prob, SSAStepper()) diff --git a/test/spatial/bracketing.jl b/test/spatial/bracketing.jl index 383e9644..31c1e23b 100644 --- a/test/spatial/bracketing.jl +++ b/test/spatial/bracketing.jl @@ -17,7 +17,7 @@ majump_rates = [0.1] # death at rate 0.1 reactstoch = [[1 => 1]] netstoch = [[1 => -1]] majump = MassActionJump(majump_rates, reactstoch, - netstoch) + netstoch) rx_rates = JP.LowHigh(JP.RxRates(n, majump)) # set up hop rates diff --git a/test/spatial/diffusion.jl b/test/spatial/diffusion.jl index 95996197..12ce485b 100644 --- a/test/spatial/diffusion.jl +++ b/test/spatial/diffusion.jl @@ -64,28 +64,28 @@ times = 0.0:(tf / num_time_points):tf algs = [NSM(), DirectCRDirect()] grids = [CartesianGridRej(dims), Graphs.grid(dims)] jump_problems = JumpProblem[JumpProblem(prob, algs[2], majumps, - hopping_constants = hopping_constants, - spatial_system = grid, - save_positions = (false, false), rng = rng) + hopping_constants = hopping_constants, + spatial_system = grid, + save_positions = (false, false), rng = rng) for grid in grids] sizehint!(jump_problems, 15 + length(jump_problems)) # flattenned jump prob push!(jump_problems, - JumpProblem(prob, NRM(), majumps, hopping_constants = hopping_constants, - spatial_system = grids[1], save_positions = (false, false), rng = rng)) + JumpProblem(prob, NRM(), majumps, hopping_constants = hopping_constants, + spatial_system = grids[1], save_positions = (false, false), rng = rng)) # hop rates of form D_s hop_constants = [hopping_rate] for alg in algs push!(jump_problems, - JumpProblem(prob, alg, majumps, hopping_constants = hop_constants, - spatial_system = grids[1], save_positions = (false, false), - rng = rng)) + JumpProblem(prob, alg, majumps, hopping_constants = hop_constants, + spatial_system = grids[1], save_positions = (false, false), + rng = rng)) push!(jump_problems, - JumpProblem(prob, alg, majumps, hopping_constants = hop_constants, - spatial_system = grids[end], save_positions = (false, false), - rng = rng)) + JumpProblem(prob, alg, majumps, hopping_constants = hop_constants, + spatial_system = grids[end], save_positions = (false, false), + rng = rng)) end # hop rates of form L_{s,i,j} @@ -93,17 +93,17 @@ hop_constants = Matrix{Vector{Float64}}(undef, size(hopping_constants)) for ci in CartesianIndices(hop_constants) (species, site) = Tuple(ci) hop_constants[ci] = repeat([hopping_constants[species, site]], - outdegree(grids[1], site)) + outdegree(grids[1], site)) end for alg in algs push!(jump_problems, - JumpProblem(prob, alg, majumps, hopping_constants = hop_constants, - spatial_system = grids[1], save_positions = (false, false), - rng = rng)) + JumpProblem(prob, alg, majumps, hopping_constants = hop_constants, + spatial_system = grids[1], save_positions = (false, false), + rng = rng)) push!(jump_problems, - JumpProblem(prob, alg, majumps, hopping_constants = hop_constants, - spatial_system = grids[end], save_positions = (false, false), - rng = rng)) + JumpProblem(prob, alg, majumps, hopping_constants = hop_constants, + spatial_system = grids[end], save_positions = (false, false), + rng = rng)) end # hop rates of form D_s * L_{i,j} @@ -114,15 +114,15 @@ for site in 1:num_nodes end for alg in algs push!(jump_problems, - JumpProblem(prob, alg, majumps, - hopping_constants = Pair(species_hop_constants, site_hop_constants), - spatial_system = grids[1], save_positions = (false, false), - rng = rng)) + JumpProblem(prob, alg, majumps, + hopping_constants = Pair(species_hop_constants, site_hop_constants), + spatial_system = grids[1], save_positions = (false, false), + rng = rng)) push!(jump_problems, - JumpProblem(prob, alg, majumps, - hopping_constants = Pair(species_hop_constants, site_hop_constants), - spatial_system = grids[end], save_positions = (false, false), - rng = rng)) + JumpProblem(prob, alg, majumps, + hopping_constants = Pair(species_hop_constants, site_hop_constants), + spatial_system = grids[end], save_positions = (false, false), + rng = rng)) end # hop rates of form D_{s,i} * L_{i,j} @@ -133,15 +133,15 @@ for site in 1:num_nodes end for alg in algs push!(jump_problems, - JumpProblem(prob, alg, majumps, - hopping_constants = Pair(species_hop_constants, site_hop_constants), - spatial_system = grids[1], save_positions = (false, false), - rng = rng)) + JumpProblem(prob, alg, majumps, + hopping_constants = Pair(species_hop_constants, site_hop_constants), + spatial_system = grids[1], save_positions = (false, false), + rng = rng)) push!(jump_problems, - JumpProblem(prob, alg, majumps, - hopping_constants = Pair(species_hop_constants, site_hop_constants), - spatial_system = grids[end], save_positions = (false, false), - rng = rng)) + JumpProblem(prob, alg, majumps, + hopping_constants = Pair(species_hop_constants, site_hop_constants), + spatial_system = grids[end], save_positions = (false, false), + rng = rng)) end # testing @@ -173,7 +173,7 @@ tspan = (0.0, 10.0) prob = DiscreteProblem(starting_state, tspan) jp = JumpProblem(prob, NSM(), majumps, hopping_constants = hopping_constants, - spatial_system = grid, save_positions = (false, false), rng = rng) + spatial_system = grid, save_positions = (false, false), rng = rng) sol = solve(jp, SSAStepper()) @test sol.u[end][1, 1] == sum(sol.u[end]) diff --git a/test/spatial/hop_rates.jl b/test/spatial/hop_rates.jl index 92a29600..56eab5a7 100644 --- a/test/spatial/hop_rates.jl +++ b/test/spatial/hop_rates.jl @@ -27,7 +27,7 @@ end normalized(distribution) = distribution / sum(distribution) function statistical_test(hop_rates, spec_propensities, target_propensities::Dict, - num_species, u, site, g, rng, rel_tol) + num_species, u, site, g, rng, rel_tol) spec_probs = normalized(spec_propensities) target_probs = normalized(target_propensities) JP.update_hop_rates!(hop_rates, 1:num_species, u, site, g) @@ -69,7 +69,7 @@ for site in 1:num_nodes target_propensities[target] = 1.0 end statistical_test(hop_rates, spec_propensities, target_propensities, num_species, u, - site, g, rng, rel_tol) + site, g, rng, rel_tol) end test_reset(hop_rates, num_nodes) @@ -85,7 +85,7 @@ for site in 1:num_nodes target_propensities[target] = 1.0 end statistical_test(hop_rates, spec_propensities, target_propensities, num_species, u, - site, g, rng, rel_tol) + site, g, rng, rel_tol) end test_reset(hop_rates, num_nodes) @@ -98,7 +98,7 @@ hop_constants[1, [1.0, 2.0], [1.0, 2.0], [1.0, 2.0, 4.0], - [1.0, 2.0], + [1.0, 2.0] ] hop_constants[2, :] = [ @@ -107,11 +107,11 @@ hop_constants[2, [3.0, 6.0], [3.0, 6.0], [3.0, 6.0, 12.0], - [3.0, 6.0], + [3.0, 6.0] ] hop_rates_structs = [ JP.HopRatesGraphDsij(hop_constants), - JP.HopRates(hop_constants, g), + JP.HopRates(hop_constants, g) ] @test hop_rates_structs[2] isa JP.HopRatesGridDsij for hop_rates in hop_rates_structs @@ -124,7 +124,7 @@ for hop_rates in hop_rates_structs for species in 1:num_species]) end statistical_test(hop_rates, spec_propensities, target_propensities, num_species, u, - site, g, rng, rel_tol) + site, g, rng, rel_tol) end end test_reset(hop_rates, num_nodes) @@ -137,11 +137,11 @@ site_hop_constants = [ [1.0, 2.0], [1.0, 2.0], [1.0, 2.0, 4.0], - [1.0, 2.0], + [1.0, 2.0] ] #[site][target_site] hop_rates_structs = [ JP.HopRatesGraphDsLij(species_hop_constants, site_hop_constants), - JP.HopRates((species_hop_constants => site_hop_constants), g), + JP.HopRates((species_hop_constants => site_hop_constants), g) ] @test hop_rates_structs[2] isa JP.HopRatesGridDsLij for hop_rates in hop_rates_structs @@ -156,7 +156,7 @@ for hop_rates in hop_rates_structs for species in 1:num_species]) end statistical_test(hop_rates, spec_propensities, target_propensities, num_species, u, - site, g, rng, rel_tol) + site, g, rng, rel_tol) end end test_reset(hop_rates, num_nodes) @@ -169,11 +169,11 @@ site_hop_constants = [ [1.0, 2.0], [1.0, 2.0], [1.0, 2.0, 4.0], - [1.0, 2.0], + [1.0, 2.0] ] #[site][target_site] hop_rates_structs = [ JP.HopRatesGraphDsiLij(species_hop_constants, site_hop_constants), - JP.HopRates((species_hop_constants => site_hop_constants), g), + JP.HopRates((species_hop_constants => site_hop_constants), g) ] @test hop_rates_structs[2] isa JP.HopRatesGridDsiLij for hop_rates in hop_rates_structs @@ -188,7 +188,7 @@ for hop_rates in hop_rates_structs for species in 1:num_species]) end statistical_test(hop_rates, spec_propensities, target_propensities, num_species, u, - site, g, rng, rel_tol) + site, g, rng, rel_tol) end end test_reset(hop_rates, num_nodes) diff --git a/test/spatial/spatial_majump.jl b/test/spatial/spatial_majump.jl index a6edc5e2..b253bad0 100644 --- a/test/spatial/spatial_majump.jl +++ b/test/spatial/spatial_majump.jl @@ -35,45 +35,47 @@ hopping_constants = [diffusivity for i in u0] # majumps uniform_majumps_1 = SpatialMassActionJump(uniform_rates[:, 1], reactstoch, netstoch) uniform_majumps_2 = SpatialMassActionJump(uniform_rates, reactstoch, netstoch) -uniform_majumps_3 = SpatialMassActionJump([1.0], reshape(uniform_rates[2, :], 1, num_nodes), - reactstoch, netstoch) # hybrid +uniform_majumps_3 = SpatialMassActionJump( + [1.0], reshape(uniform_rates[2, :], 1, num_nodes), + reactstoch, netstoch) # hybrid uniform_majumps_4 = SpatialMassActionJump(MassActionJump(uniform_rates[:, 1], reactstoch, - netstoch)) + netstoch)) uniform_majumps = [ uniform_majumps_1, uniform_majumps_2, uniform_majumps_3, - uniform_majumps_4, + uniform_majumps_4 ] non_uniform_majumps_1 = SpatialMassActionJump(non_uniform_rates, reactstoch, netstoch) # reactions are zero outside of center site non_uniform_majumps_2 = SpatialMassActionJump([1.0], - reshape(non_uniform_rates[2, :], 1, - num_nodes), reactstoch, netstoch) # birth everywhere, death only at center site -non_uniform_majumps_3 = SpatialMassActionJump([1.0 0.0 0.0 0.0 0.0; - 0.0 0.0 0.0 0.0 death_rate], reactstoch, - netstoch) # birth on the left, death on the right + reshape(non_uniform_rates[2, :], 1, + num_nodes), reactstoch, netstoch) # birth everywhere, death only at center site +non_uniform_majumps_3 = SpatialMassActionJump( + [1.0 0.0 0.0 0.0 0.0; + 0.0 0.0 0.0 0.0 death_rate], reactstoch, + netstoch) # birth on the left, death on the right non_uniform_majumps = [non_uniform_majumps_1, non_uniform_majumps_2, non_uniform_majumps_3] # put together the JumpProblem's uniform_jump_problems = JumpProblem[JumpProblem(prob, NSM(), majump, - hopping_constants = hopping_constants, - spatial_system = grid, - save_positions = (false, false), rng = rng) + hopping_constants = hopping_constants, + spatial_system = grid, + save_positions = (false, false), rng = rng) for majump in uniform_majumps] # flattenned append!(uniform_jump_problems, - JumpProblem[JumpProblem(prob, NRM(), majump, hopping_constants = hopping_constants, - spatial_system = grid, save_positions = (false, false), - rng = rng) - for majump in uniform_majumps]) + JumpProblem[JumpProblem(prob, NRM(), majump, hopping_constants = hopping_constants, + spatial_system = grid, save_positions = (false, false), + rng = rng) + for majump in uniform_majumps]) # non-uniform non_uniform_jump_problems = JumpProblem[JumpProblem(prob, NSM(), majump, - hopping_constants = hopping_constants, - spatial_system = grid, - save_positions = (false, false), - rng = rng) + hopping_constants = hopping_constants, + spatial_system = grid, + save_positions = (false, false), + rng = rng) for majump in non_uniform_majumps] # testing diff --git a/test/spatial/topology.jl b/test/spatial/topology.jl index c1bd7519..9c93445b 100644 --- a/test/spatial/topology.jl +++ b/test/spatial/topology.jl @@ -17,7 +17,7 @@ num_samples = 10^5 rel_tol = 0.01 grids = [ JP.CartesianGridRej(dims), - Graphs.grid(dims), + Graphs.grid(dims) ] for grid in grids show(io, "text/plain", grid) diff --git a/test/splitcoupled.jl b/test/splitcoupled.jl index ff22adba..e871e2c0 100644 --- a/test/splitcoupled.jl +++ b/test/splitcoupled.jl @@ -16,8 +16,9 @@ jump_prob = JumpProblem(prob, Direct(), jump1; rng = rng) jump_prob_control = JumpProblem(prob_control, Direct(), jump1; rng = rng) coupling_map = [(1, 1)] -coupled_prob = SplitCoupledJumpProblem(jump_prob, jump_prob_control, Direct(), coupling_map; - rng = rng) +coupled_prob = SplitCoupledJumpProblem( + jump_prob, jump_prob_control, Direct(), coupling_map; + rng = rng) @time sol = solve(coupled_prob, FunctionMap()) @time solve(jump_prob, FunctionMap()) @@ -43,8 +44,9 @@ prob = ODEProblem(f, [1.0], (0.0, 1.0)) prob_control = ODEProblem(f, [1.0], (0.0, 1.0)) jump_prob = JumpProblem(prob, Direct(), jump1; rng = rng) jump_prob_control = JumpProblem(prob_control, Direct(), jump2; rng = rng) -coupled_prob = SplitCoupledJumpProblem(jump_prob, jump_prob_control, Direct(), coupling_map; - rng = rng) +coupled_prob = SplitCoupledJumpProblem( + jump_prob, jump_prob_control, Direct(), coupling_map; + rng = rng) sol = solve(coupled_prob, Tsit5()) @test mean([abs(s[1] - s[2]) for s in sol.u]) <= 5.0 @@ -53,8 +55,9 @@ prob = SDEProblem(f, g, [1.0], (0.0, 1.0)) prob_control = SDEProblem(f, g, [1.0], (0.0, 1.0)) jump_prob = JumpProblem(prob, Direct(), jump1; rng = rng) jump_prob_control = JumpProblem(prob_control, Direct(), jump1; rng = rng) -coupled_prob = SplitCoupledJumpProblem(jump_prob, jump_prob_control, Direct(), coupling_map; - rng = rng) +coupled_prob = SplitCoupledJumpProblem( + jump_prob, jump_prob_control, Direct(), coupling_map; + rng = rng) sol = solve(coupled_prob, SRIW1()) @test mean([abs(s[1] - s[2]) for s in sol.u]) <= 5.0 @@ -63,8 +66,9 @@ prob = ODEProblem(f, [1.0], (0.0, 1.0)) prob_control = SDEProblem(f, g, [1.0], (0.0, 1.0)) jump_prob = JumpProblem(prob, Direct(), jump1; rng = rng) jump_prob_control = JumpProblem(prob_control, Direct(), jump1; rng = rng) -coupled_prob = SplitCoupledJumpProblem(jump_prob, jump_prob_control, Direct(), coupling_map; - rng = rng) +coupled_prob = SplitCoupledJumpProblem( + jump_prob, jump_prob_control, Direct(), coupling_map; + rng = rng) sol = solve(coupled_prob, SRIW1()) @test mean([abs(s[1] - s[2]) for s in sol.u]) <= 5.0 @@ -77,8 +81,9 @@ prob = DiscreteProblem([1.0], (0.0, 1.0)) prob_control = SDEProblem(f, g, [1.0], (0.0, 1.0)) jump_prob = JumpProblem(prob, Direct(), jump1; rng = rng) jump_prob_control = JumpProblem(prob_control, Direct(), jump1; rng = rng) -coupled_prob = SplitCoupledJumpProblem(jump_prob, jump_prob_control, Direct(), coupling_map; - rng = rng) +coupled_prob = SplitCoupledJumpProblem( + jump_prob, jump_prob_control, Direct(), coupling_map; + rng = rng) sol = solve(coupled_prob, SRIW1()) # test mass action jumps coupled to ODE @@ -92,7 +97,7 @@ f = function (du, u, p, t) end odeprob = ODEProblem(f, [10.0], (0.0, 10.0)) jump_prob = JumpProblem(odeprob, Direct(), majumps, save_positions = (false, false); - rng = rng) + rng = rng) Nsims = 8000 Amean = 0.0 for i in 1:Nsims diff --git a/test/ssa_callback_test.jl b/test/ssa_callback_test.jl index d6b004d5..2fb0e23f 100644 --- a/test/ssa_callback_test.jl +++ b/test/ssa_callback_test.jl @@ -48,7 +48,7 @@ finalizer_called = 0 fuel_finalize(cb, u, t, integrator) = global finalizer_called += 1 cb2 = DiscreteCallback(condition, fuel_affect!, initialize = fuel_init!, - finalize = fuel_finalize) + finalize = fuel_finalize) sol = solve(jump_prob, SSAStepper(), callback = cb2) for tstop in random_tstops @test tstop ∈ sol.t @@ -71,7 +71,7 @@ function paffect!(integrator) reset_aggregated_jumps!(integrator) end sol = solve(jprob, SSAStepper(), tstops = [1000.0], - callback = DiscreteCallback(pcondit, paffect!)) + callback = DiscreteCallback(pcondit, paffect!)) @test all(p .== [0.0, 1.0]) @test sol[1, end] == 100 @@ -80,7 +80,7 @@ maj1 = MassActionJump([1 => 1], [1 => -1, 2 => 1]; param_idxs = 1) maj2 = MassActionJump([2 => 1], [1 => 1, 2 => -1]; param_idxs = 2) jprob = JumpProblem(dprob, Direct(), maj1, maj2, save_positions = (false, false), rng = rng) sol = solve(jprob, SSAStepper(), tstops = [1000.0], - callback = DiscreteCallback(pcondit, paffect!)) + callback = DiscreteCallback(pcondit, paffect!)) @test all(p .== [0.0, 1.0]) @test sol[1, end] == 100 @@ -88,27 +88,27 @@ p2 = [1.0, 0.0, 0.0] maj3 = MassActionJump([1 => 1], [1 => -1, 2 => 1]; param_idxs = 3) dprob = DiscreteProblem(u₀, tspan, p2) jprob = JumpProblem(dprob, Direct(), maj1, maj2, maj3, save_positions = (false, false), - rng = rng) + rng = rng) sol = solve(jprob, SSAStepper(), tstops = [1000.0], - callback = DiscreteCallback(pcondit, paffect!)) + callback = DiscreteCallback(pcondit, paffect!)) @test all(p2 .== [0.0, 1.0, 0.0]) @test sol[1, end] == 100 p2 .= [1.0, 0.0, 0.0] jprob = JumpProblem(dprob, Direct(), JumpSet(; massaction_jumps = [maj1, maj2, maj3]), - save_positions = (false, false), rng = rng) + save_positions = (false, false), rng = rng) sol = solve(jprob, SSAStepper(), tstops = [1000.0], - callback = DiscreteCallback(pcondit, paffect!)) + callback = DiscreteCallback(pcondit, paffect!)) @test all(p2 .== [0.0, 1.0, 0.0]) @test sol[1, end] == 100 p .= [1.0, 0.0] dprob = DiscreteProblem(u₀, tspan, p) maj4 = MassActionJump([[1 => 1], [2 => 1]], [[1 => -1, 2 => 1], [1 => 1, 2 => -1]]; - param_idxs = [1, 2]) + param_idxs = [1, 2]) jprob = JumpProblem(dprob, Direct(), maj4, save_positions = (false, false), rng = rng) sol = solve(jprob, SSAStepper(), tstops = [1000.0], - callback = DiscreteCallback(pcondit, paffect!)) + callback = DiscreteCallback(pcondit, paffect!)) @test all(p .== [0.0, 1.0]) @test sol[1, end] == 100 @@ -119,12 +119,12 @@ maj5 = MassActionJump([[1 => 2]], [[1 => -1, 2 => 1]]; param_idxs = [1]) jprob = JumpProblem(dprob, Direct(), maj5, save_positions = (false, false), rng = rng) @test all(jprob.massaction_jump.scaled_rates .== [0.5]) jprob = JumpProblem(dprob, Direct(), maj5, save_positions = (false, false), rng = rng, - scale_rates = false) + scale_rates = false) @test all(jprob.massaction_jump.scaled_rates .== [1.0]) # test for https://github.com/SciML/JumpProcesses.jl/issues/239 maj6 = MassActionJump([[1 => 1], [2 => 1]], [[1 => -1, 2 => 1], [1 => 1, 2 => -1]]; - param_idxs = [1, 2]) + param_idxs = [1, 2]) p = (0.1, 0.1) dprob = DiscreteProblem([10, 0], (0.0, 100.0), p) jprob = JumpProblem(dprob, Direct(), maj6; save_positions = (false, false), rng = rng) diff --git a/test/ssa_tests.jl b/test/ssa_tests.jl index 66fc14b8..e82c5095 100644 --- a/test/ssa_tests.jl +++ b/test/ssa_tests.jl @@ -38,7 +38,7 @@ sol = solve(jump_prob, SSAStepper(), save_start = false) @test sol.t[end] == 3.0 jump_prob = JumpProblem(prob, Direct(), jump, jump2, save_positions = (false, false); - rng = rng) + rng = rng) sol = solve(jump_prob, SSAStepper(), save_start = false, save_end = false) @test isempty(sol.t) && isempty(sol.u) diff --git a/test/thread_safety.jl b/test/thread_safety.jl index d78430ae..4f1db556 100644 --- a/test/thread_safety.jl +++ b/test/thread_safety.jl @@ -12,4 +12,4 @@ dprob = DiscreteProblem(u0, tspan, params) jprob = JumpProblem(dprob, Direct(), maj; rng = rng) solve(EnsembleProblem(jprob), SSAStepper(), EnsembleThreads(); trajectories = 10) solve(EnsembleProblem(jprob; safetycopy = true), SSAStepper(), EnsembleThreads(); - trajectories = 10) + trajectories = 10) diff --git a/test/variable_rate.jl b/test/variable_rate.jl index 232ff067..83de0f46 100644 --- a/test/variable_rate.jl +++ b/test/variable_rate.jl @@ -125,9 +125,9 @@ prob = ODEProblem(f3, [1.0 2.0; 3.0 4.0], (0.0, 1.0)) rate3(u, p, t) = u[1] + u[2] function affect3!(integrator) (integrator.u[1] = 0.25; - integrator.u[2] = 0.5; - integrator.u[3] = 0.75; - integrator.u[4] = 1) + integrator.u[2] = 0.5; + integrator.u[3] = 0.75; + integrator.u[4] = 1) end jump = VariableRateJump(rate3, affect3!) jump_prob = JumpProblem(prob, Direct(), jump; rng = rng) @@ -173,7 +173,7 @@ let react_stoich_ = [Vector{Pair{Int, Int}}()] net_stoich_ = [[1 => 1]] mass_action_jump_ = MassActionJump(maj_rate, react_stoich_, net_stoich_; - scale_rates = false) + scale_rates = false) affect! = function (integrator) integrator.u[1] -= 1 @@ -187,12 +187,12 @@ let tspan = (0.0, 30.0) dprob_ = DiscreteProblem(u0, tspan) @test_throws ErrorException JumpProblem(dprob_, alg, jumpset_, - save_positions = (false, false)) + save_positions = (false, false)) vrj = VariableRateJump(cs_rate1, affect!; urate = ((u, p, t) -> 1.0), - rateinterval = ((u, p, t) -> 1.0)) + rateinterval = ((u, p, t) -> 1.0)) @test_throws ErrorException JumpProblem(dprob_, alg, mass_action_jump_, vrj; - save_positions = (false, false)) + save_positions = (false, false)) end end @@ -223,7 +223,7 @@ let end test_jump = VariableRateJump(test_rate, test_affect!; urate = test_urate, - rateinterval = (u, p, t) -> 1.0) + rateinterval = (u, p, t) -> 1.0) dprob = DiscreteProblem([0], (0.0, 1.0), nothing) jprob = JumpProblem(dprob, Coevolve(), test_jump; dep_graph = [[1]])