Publish refreshed Jumps benchmarks

script/Jumps/Diffusion_CTRW.jl

@@ -1,29 +1,33 @@
-
 using Catalyst, JumpProcesses, JumpProblemLibrary, Plots, Statistics, DataFrames
 
 
 N = 256
 h = 1 / N
-u0 = 10 * ones(Int64, N)
-tf = .01
-methods = (Direct(), FRM(), SortingDirect(), NRM(), DirectCR(), RSSA(), RSSACR(), Coevolve(), RDirect())
-shortlabels = [string(leg)[15:end-2] for leg in methods]
+tf = 0.01
+methods = (
+    Direct(), FRM(), SortingDirect(), NRM(), DirectCR(),
+    RSSA(), RSSACR(), Coevolve(), RDirect(),
+)
+shortlabels = [string(leg)[15:(end - 2)] for leg in methods]
 jprob = JumpProblemLibrary.prob_jump_diffnetwork
 rn = jprob.network(N)
-prob = DiscreteProblem(rn, u0, (0.0, tf), [1 / (h*h)])
-ploth = plot(reuse=false)
-for (i,method) in enumerate(methods)
+u0 = jprob.u0(rn, N)
+rates = [:K => 1 / (h * h)]
+ploth = plot(reuse = false)
+for (i, method) in enumerate(methods)
     println("Benchmarking method: ", method)
-    jump_prob = JumpProblem(rn, prob, method, save_positions=(false, false))
-    sol = solve(jump_prob, SSAStepper(); saveat=tf/1000.)
-    plot!(ploth, sol.t, sol[Int(N//2),:], label=shortlabels[i])
+    jump_prob = JumpProblem(
+        rn, u0, (0.0, tf), rates; aggregator = method, save_positions = (false, false)
+    )
+    sol = solve(jump_prob, SSAStepper(); saveat = tf / 1000.0)
+    plot!(ploth, sol.t, sol[Int(N // 2), :], label = shortlabels[i])
 end
-plot!(ploth, title="Population at middle lattice site", xlabel="time")
+plot!(ploth, title = "Population at middle lattice site", xlabel = "time")
 
 
 function run_benchmark!(t, jump_prob, stepper)
     sol = solve(jump_prob, stepper)
-    @inbounds for i in 1:length(t)
+    return @inbounds for i in 1:length(t)
         t[i] = @elapsed (sol = solve(jump_prob, stepper))
     end
 end
@@ -32,34 +36,37 @@ end
 nsims = 50
 benchmarks = Vector{Vector{Float64}}()
 for method in methods
-    jump_prob = JumpProblem(rn, prob, method, save_positions=(false, false))
+    jump_prob = JumpProblem(
+        rn, u0, (0.0, tf), rates; aggregator = method, save_positions = (false, false)
+    )
     stepper = SSAStepper()
     t = Vector{Float64}(undef, nsims)
     run_benchmark!(t, jump_prob, stepper)
     push!(benchmarks, t)
 end
 
 
-medtimes = Vector{Float64}(undef,length(methods))
-stdtimes = Vector{Float64}(undef,length(methods))
-avgtimes = Vector{Float64}(undef,length(methods))
+medtimes = Vector{Float64}(undef, length(methods))
+stdtimes = Vector{Float64}(undef, length(methods))
+avgtimes = Vector{Float64}(undef, length(methods))
 for i in 1:length(methods)
     medtimes[i] = median(benchmarks[i])
     avgtimes[i] = mean(benchmarks[i])
     stdtimes[i] = std(benchmarks[i])
 end
 
-df = DataFrame(names=shortlabels, medtimes=medtimes, relmedtimes=(medtimes/medtimes[1]),
-               avgtimes=avgtimes, std=stdtimes, cv=stdtimes./avgtimes)
+df = DataFrame(
+    names = shortlabels, medtimes = medtimes, relmedtimes = (medtimes / medtimes[1]),
+    avgtimes = avgtimes, std = stdtimes, cv = stdtimes ./ avgtimes
+)
 
 
-sa = [string(round(mt,digits=4),"s") for mt in df.medtimes]
-bar(df.names, df.relmedtimes, legend=:false)
-scatter!(df.names, .05 .+ df.relmedtimes, markeralpha=0, series_annotations=sa)
+sa = [string(round(mt, digits = 4), "s") for mt in df.medtimes]
+bar(df.names, df.relmedtimes, legend = :false)
+scatter!(df.names, 0.05 .+ df.relmedtimes, markeralpha = 0, series_annotations = sa)
 ylabel!("median relative to Direct")
 title!("256 Site 1D Diffusion CTRW")
 
 
 using SciMLBenchmarks
-SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])
-
+SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder], WEAVE_ARGS[:file])

script/Jumps/EGFR_Benchmark.jl

@@ -1,45 +1,54 @@
+using JumpProcesses, Plots, StableRNGs, BenchmarkTools, ReactionNetworkImporters,
+    StatsPlots, Catalyst
 
-using JumpProcesses, Plots, StableRNGs, BenchmarkTools, ReactionNetworkImporters, StatsPlots, Catalyst
 
-
-tf = 12.
+tf = 12.0
 rng = StableRNG(53124)
 algs = [NRM(), CCNRM(), DirectCR(), RSSACR()]
-egfr_net = loadrxnetwork(BNGNetwork(), joinpath(@__DIR__, "Data/egfr_net.net"));
-dprob = DiscreteProblem(complete(egfr_net.rn), egfr_net.u0, (0., tf), egfr_net.p)
-dprob = remake(dprob,u0=Int64.(dprob.u0))
+egfr_net = complete(loadrxnetwork(BNGNetwork(), joinpath(@__DIR__, "Data/egfr_net.net")))
+egfr_u0 = Catalyst.get_u0_map(egfr_net)
+egfr_parameters = Catalyst.get_parameter_map(egfr_net)
 
-plt = plot(title="Dimer concentrations")
+plt = plot(title = "Dimer concentrations")
 for alg in algs
-    jprob = JumpProblem(complete(egfr_net.rn), dprob, alg)
-    sol = solve(jprob, SSAStepper(), saveat = tf/200)
-    plot!(plt, sol, idxs = :Dimers, label="$alg")
+    jprob = JumpProblem(
+        egfr_net, egfr_u0, (0.0, tf), egfr_parameters;
+        aggregator = alg, u0_eltype = Int64
+    )
+    sol = solve(jprob, SSAStepper(), saveat = tf / 200)
+    plot!(plt, sol, idxs = :Dimers, label = "$alg")
 end
 plot!(plt)
 
 
 function benchmark_and_bar_plot(model, end_time, algs)
     times = Vector{Float64}()
-    alg_names = ["$s"[15:end-2] for s in algs]
+    alg_names = ["$s"[15:(end - 2)] for s in algs]
+    u0 = Catalyst.get_u0_map(model)
+    parameters = Catalyst.get_parameter_map(model)
 
     benchmarks = Vector{BenchmarkTools.Trial}(undef, length(algs))
     for (i, alg) in enumerate(algs)
         alg_name = alg_names[i]
         println("Benchmarking $alg_name")
-        dprob = DiscreteProblem(complete(model.rn), model.u0, (0., end_time), model.p)
-        dprob = remake(dprob,u0 = Int64.(dprob.u0))
-        jprob = JumpProblem(complete(model.rn), dprob, alg; rng, save_positions = (false, false))
+        jprob = JumpProblem(
+            model, u0, (0.0, end_time), parameters;
+            aggregator = alg, rng, save_positions = (false, false), u0_eltype = Int64
+        )
 
         b = @benchmarkable solve($jprob; saveat = $end_time) samples = 5 seconds = 7200
         bm = run(b)
-        push!(times, median(bm).time/1e9)
+        push!(times, median(bm).time / 1.0e9)
     end
 
-    bar(alg_names, times, xlabel = "Algorithm", ylabel = "Average Time (s)", title = "SSA Runtime for EGFR network", legend = false)
+    return bar(
+        alg_names, times, xlabel = "Algorithm", ylabel = "Average Time (s)",
+        title = "SSA Runtime for EGFR network", legend = false
+    )
 end
 
 
-tf = 12. 
+tf = 12.0
 rng = StableRNG(53124)
 algs = [NRM(), CCNRM(), DirectCR(), RSSACR()]
 
@@ -49,4 +58,3 @@ plt
 
 using SciMLBenchmarks
 SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder], WEAVE_ARGS[:file])
-

script/Jumps/Mendes_multistate_example.jl

@@ -1,52 +1,56 @@
-
 using Catalyst, JumpProcesses, JumpProblemLibrary, Plots, Statistics
 fmt = :png
 
 
 jprob = JumpProblemLibrary.prob_jump_multistate
-rn  = jprob.network
+rn = jprob.network
 reactions(rn)
 
 
-methods = (Direct(), FRM(), SortingDirect(), NRM(), DirectCR(), RSSA(), RSSACR(), Coevolve(), RDirect())
-shortlabels = [string(leg)[15:end-2] for leg in methods]
-tf      = 10.0 * jprob.tstop
-prob    = DiscreteProblem(rn, jprob.u0, (0.0, tf), jprob.rates)
+methods = (
+    Direct(), FRM(), SortingDirect(), NRM(), DirectCR(),
+    RSSA(), RSSACR(), Coevolve(), RDirect(),
+)
+shortlabels = [string(leg)[15:(end - 2)] for leg in methods]
+tf = 10.0 * jprob.tstop
 varlegs = ["A_P" "A_bound_P" "A_unbound_P" "RLA_P"]
 @variables t S7(t) S8(t) S9(t)
 varsyms = [
-    [S7,S8,S9],
+    [S7, S8, S9],
     [S9],
-    [S7,S8],
-    [S7]
+    [S7, S8],
+    [S7],
 ]
 varidxs = []
 for vars in varsyms
-    push!(varidxs, [findfirst(isequal(sym),unknowns(rn)) for sym in vars])
+    push!(varidxs, [findfirst(isequal(sym), unknowns(rn)) for sym in vars])
 end
 
 
 p = []
-for (i,method) in enumerate(methods)
-    jump_prob = JumpProblem(rn, prob, method, save_positions=(false, false))
-    sol = solve(jump_prob, SSAStepper(), saveat=tf/1000.)
+for (i, method) in enumerate(methods)
+    jump_prob = JumpProblem(
+        rn, jprob.u0, (0.0, tf), jprob.rates;
+        aggregator = method, save_positions = (false, false)
+    )
+    sol = solve(jump_prob, SSAStepper(), saveat = tf / 1000.0)
     solv = zeros(1001, 4)
-    for (i,varidx) in enumerate(varidxs)
-        solv[:,i] = sum(sol[varidx,:], dims=1)
+    for (i, varidx) in enumerate(varidxs)
+        solv[:, i] = sum(sol[varidx, :], dims = 1)
     end
     if i < length(methods)
-        push!(p, plot(sol.t, solv, title=shortlabels[i], legend=false, format=fmt))
+        push!(p, plot(sol.t, solv, title = shortlabels[i], legend = false, format = fmt))
     else
-        push!(p, plot(sol.t, solv, title=shortlabels[i], legend=false, format=fmt))
+        push!(p, plot(sol.t, solv, title = shortlabels[i], legend = false, format = fmt))
     end
 end
-push!(p, plot((1:4)', framestyle = :none, legend=:inside, labels=varlegs))
-plot(p..., layout=(6,2), format=fmt)
+push!(p, plot((1:4)', framestyle = :none, legend = :inside, labels = varlegs))
+plot(p..., layout = (6, 2), format = fmt)
 
 
 function run_benchmark!(t, jump_prob, stepper)
     sol = solve(jump_prob, stepper)
-    @inbounds for i in 1:length(t)
+    return @inbounds for i in 1:length(t)
         t[i] = @elapsed (sol = solve(jump_prob, stepper))
     end
 end
@@ -55,34 +59,40 @@ end
 nsims = 100
 benchmarks = Vector{Vector{Float64}}()
 for method in methods
-    jump_prob = JumpProblem(rn, prob, method, save_positions=(false, false))
+    jump_prob = JumpProblem(
+        rn, jprob.u0, (0.0, tf), jprob.rates;
+        aggregator = method, save_positions = (false, false)
+    )
     stepper = SSAStepper()
     time = Vector{Float64}(undef, nsims)
     run_benchmark!(time, jump_prob, stepper)
     push!(benchmarks, time)
 end
 
 
-medtimes = Vector{Float64}(undef,length(methods))
-stdtimes = Vector{Float64}(undef,length(methods))
-avgtimes = Vector{Float64}(undef,length(methods))
+medtimes = Vector{Float64}(undef, length(methods))
+stdtimes = Vector{Float64}(undef, length(methods))
+avgtimes = Vector{Float64}(undef, length(methods))
 for i in 1:length(methods)
     medtimes[i] = median(benchmarks[i])
     avgtimes[i] = mean(benchmarks[i])
     stdtimes[i] = std(benchmarks[i])
 end
 using DataFrames
 
-df = DataFrame(names=shortlabels, medtimes=medtimes, relmedtimes=(medtimes/medtimes[1]),
-               avgtimes=avgtimes, std=stdtimes, cv=stdtimes./avgtimes)
+df = DataFrame(
+    names = shortlabels, medtimes = medtimes, relmedtimes = (medtimes / medtimes[1]),
+    avgtimes = avgtimes, std = stdtimes, cv = stdtimes ./ avgtimes
+)
 
-sa = [text(string(round(mt,digits=3),"s"),:center,12) for mt in df.medtimes]
-bar(df.names,df.relmedtimes,legend=:false, fmt=fmt)
-scatter!(df.names, .05 .+ df.relmedtimes, markeralpha=0, series_annotations=sa, fmt=fmt)
+sa = [text(string(round(mt, digits = 3), "s"), :center, 12) for mt in df.medtimes]
+bar(df.names, df.relmedtimes, legend = :false, fmt = fmt)
+scatter!(
+    df.names, 0.05 .+ df.relmedtimes, markeralpha = 0, series_annotations = sa, fmt = fmt
+)
 ylabel!("median relative to Direct")
 title!("Multistate Model")
 
 
 using SciMLBenchmarks
-SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])
-
+SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder], WEAVE_ARGS[:file])