diff --git a/benchmarks/HybridJumps/Manifest.toml b/benchmarks/HybridJumps/Manifest.toml index da875a8bf..6bda27310 100644 --- a/benchmarks/HybridJumps/Manifest.toml +++ b/benchmarks/HybridJumps/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.10.9" manifest_format = "2.0" -project_hash = "59bb9ca342edc40032e032bc047178bab9947f84" +project_hash = "655ee19d9fa65ef20fe7edfa85c2bffd1ce645f6" [[deps.ADTypes]] git-tree-sha1 = "e2478490447631aedba0823d4d7a80b2cc8cdb32" @@ -548,6 +548,17 @@ weakdeps = ["ChainRulesCore", "EnzymeCore"] DispatchDoctorChainRulesCoreExt = "ChainRulesCore" DispatchDoctorEnzymeCoreExt = "EnzymeCore" +[[deps.Distances]] +deps = ["LinearAlgebra", "Statistics", "StatsAPI"] +git-tree-sha1 = "c7e3a542b999843086e2f29dac96a618c105be1d" +uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" +version = "0.10.12" +weakdeps = ["ChainRulesCore", "SparseArrays"] + + [deps.Distances.extensions] + DistancesChainRulesCoreExt = "ChainRulesCore" + DistancesSparseArraysExt = "SparseArrays" + [[deps.Distributed]] deps = ["Random", "Serialization", "Sockets"] uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" @@ -569,9 +580,9 @@ version = "0.25.120" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [[deps.DocStringExtensions]] -git-tree-sha1 = "e7b7e6f178525d17c720ab9c081e4ef04429f860" +git-tree-sha1 = "7442a5dfe1ebb773c29cc2962a8980f47221d76c" uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -version = "0.9.4" +version = "0.9.5" [[deps.DomainSets]] deps = ["CompositeTypes", "IntervalSets", "LinearAlgebra", "Random", "StaticArrays"] @@ -620,9 +631,9 @@ uuid = "4e289a0a-7415-4d19-859d-a7e5c4648b56" version = "1.0.5" [[deps.EnzymeCore]] -git-tree-sha1 = "1eb59f40a772d0fbd4cb75e00b3fa7f5f79c975a" +git-tree-sha1 = "7d7822a643c33bbff4eab9c87ca8459d7c688db0" uuid = "f151be2c-9106-41f4-ab19-57ee4f262869" -version = "0.8.9" +version = "0.8.11" weakdeps = ["Adapt"] [deps.EnzymeCore.extensions] @@ -908,10 +919,10 @@ uuid = "3b182d85-2403-5c21-9c21-1e1f0cc25472" version = "1.3.15+0" [[deps.Graphs]] -deps = ["ArnoldiMethod", "Compat", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"] -git-tree-sha1 = "3169fd3440a02f35e549728b0890904cfd4ae58a" +deps = ["ArnoldiMethod", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"] +git-tree-sha1 = "c5abfa0ae0aaee162a3fbb053c13ecda39be545b" uuid = "86223c79-3864-5bf0-83f7-82e725a168b6" -version = "1.12.1" +version = "1.13.0" [[deps.Grisu]] git-tree-sha1 = "53bb909d1151e57e2484c3d1b53e19552b887fb2" @@ -944,9 +955,9 @@ version = "0.3.28" [[deps.IJulia]] deps = ["Base64", "Conda", "Dates", "InteractiveUtils", "JSON", "Logging", "Markdown", "MbedTLS", "Pkg", "Printf", "REPL", "Random", "SoftGlobalScope", "UUIDs", "ZMQ"] -git-tree-sha1 = "be30be76e25b0aff2c9a85930ed3ac34c5f10c83" +git-tree-sha1 = "5bea4d841d5db750a1bb5ccdc37b2759dfa80066" uuid = "7073ff75-c697-5162-941a-fcdaad2a7d2a" -version = "1.27.0" +version = "1.28.1" [[deps.IfElse]] git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1" @@ -1077,9 +1088,9 @@ version = "0.4.10" [[deps.JumpProcesses]] deps = ["ArrayInterface", "DataStructures", "DiffEqBase", "DiffEqCallbacks", "DocStringExtensions", "FunctionWrappers", "Graphs", "LinearAlgebra", "Markdown", "PoissonRandom", "Random", "RandomNumbers", "RecursiveArrayTools", "Reexport", "SciMLBase", "Setfield", "StaticArrays", "SymbolicIndexingInterface", "UnPack"] -git-tree-sha1 = "216c196df09c8b80a40a2befcb95760eb979bcfd" +git-tree-sha1 = "fb7fd516de38db80f50fe15e57d44da2836365e7" uuid = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5" -version = "9.15.0" +version = "9.16.0" weakdeps = ["FastBroadcast"] [[deps.Krylov]] @@ -1181,6 +1192,12 @@ git-tree-sha1 = "4c4a5a0bdac384fd0e2fc9969ad600e4473868be" uuid = "b4f0291d-fe17-52bc-9479-3d1a343d9043" version = "4.0.1" +[[deps.LevyArea]] +deps = ["LinearAlgebra", "Random", "SpecialFunctions"] +git-tree-sha1 = "56513a09b8e0ae6485f34401ea9e2f31357958ec" +uuid = "2d8b4e74-eb68-11e8-0fb9-d5eb67b50637" +version = "1.0.0" + [[deps.LibCURL]] deps = ["LibCURL_jll", "MozillaCACerts_jll"] uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21" @@ -1266,9 +1283,9 @@ uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" [[deps.LinearSolve]] deps = ["ArrayInterface", "ChainRulesCore", "ConcreteStructs", "DocStringExtensions", "EnumX", "GPUArraysCore", "InteractiveUtils", "Krylov", "LazyArrays", "Libdl", "LinearAlgebra", "MKL_jll", "Markdown", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "Setfield", "StaticArraysCore", "UnPack"] -git-tree-sha1 = "c618a6a774d5712c6bf02dbcceb51b6dc6b9bb89" +git-tree-sha1 = "c0d1a91a50af6778863d320761f807f641f74935" uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" -version = "3.16.0" +version = "3.17.0" [deps.LinearSolve.extensions] LinearSolveBandedMatricesExt = "BandedMatrices" @@ -1358,9 +1375,9 @@ uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" [[deps.MathOptInterface]] deps = ["BenchmarkTools", "CodecBzip2", "CodecZlib", "DataStructures", "ForwardDiff", "JSON3", "LinearAlgebra", "MutableArithmetics", "NaNMath", "OrderedCollections", "PrecompileTools", "Printf", "SparseArrays", "SpecialFunctions", "Test"] -git-tree-sha1 = "d4eb6037ce33f974a58a2ea9ccba28beab2a93c1" +git-tree-sha1 = "2f2c18c6acab9042557bdb0af8c3a14dd7b64413" uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" -version = "1.40.1" +version = "1.41.0" [[deps.MaybeInplace]] deps = ["ArrayInterface", "LinearAlgebra", "MacroTools"] @@ -1399,9 +1416,9 @@ uuid = "a63ad114-7e13-5084-954f-fe012c677804" [[deps.ModelingToolkit]] deps = ["ADTypes", "AbstractTrees", "ArrayInterface", "BlockArrays", "ChainRulesCore", "Combinatorics", "CommonSolve", "Compat", "ConstructionBase", "DataStructures", "DiffEqBase", "DiffEqCallbacks", "DiffEqNoiseProcess", "DiffRules", "DifferentiationInterface", "Distributed", "Distributions", "DocStringExtensions", "DomainSets", "DynamicQuantities", "EnumX", "ExprTools", "FindFirstFunctions", "ForwardDiff", "FunctionWrappers", "FunctionWrappersWrappers", "Graphs", "InteractiveUtils", "JuliaFormatter", "JumpProcesses", "Latexify", "Libdl", "LinearAlgebra", "MLStyle", "Moshi", "NaNMath", "NonlinearSolve", "OffsetArrays", "OrderedCollections", "PrecompileTools", "RecursiveArrayTools", "Reexport", "RuntimeGeneratedFunctions", "SCCNonlinearSolve", "SciMLBase", "SciMLStructures", "Serialization", "Setfield", "SimpleNonlinearSolve", "SparseArrays", "SpecialFunctions", "StaticArrays", "SymbolicIndexingInterface", "SymbolicUtils", "Symbolics", "URIs", "UnPack", "Unitful"] -git-tree-sha1 = "7be3c50f3c384293aad704626c1feeca3d295ed7" +git-tree-sha1 = "efbc9ec328e0b3c0353bc8afc0a811e7d616bade" uuid = "961ee093-0014-501f-94e3-6117800e7a78" -version = "9.80.2" +version = "9.80.5" [deps.ModelingToolkit.extensions] MTKBifurcationKitExt = "BifurcationKit" @@ -1442,9 +1459,9 @@ version = "0.5.9" [[deps.Mustache]] deps = ["Printf", "Tables"] -git-tree-sha1 = "3b2db451a872b20519ebb0cec759d3d81a1c6bcb" +git-tree-sha1 = "3cbd5dda543bc59f2e482607ccf84b633724fc32" uuid = "ffc61752-8dc7-55ee-8c37-f3e9cdd09e70" -version = "1.0.20" +version = "1.0.21" [[deps.MutableArithmetics]] deps = ["LinearAlgebra", "SparseArrays", "Test"] @@ -1454,9 +1471,15 @@ version = "1.6.4" [[deps.NLSolversBase]] deps = ["ADTypes", "DifferentiationInterface", "Distributed", "FiniteDiff", "ForwardDiff"] -git-tree-sha1 = "b14c7be6046e7d48e9063a0053f95ee0fc954176" +git-tree-sha1 = "25a6638571a902ecfb1ae2a18fc1575f86b1d4df" uuid = "d41bc354-129a-5804-8e4c-c37616107c6c" -version = "7.9.1" +version = "7.10.0" + +[[deps.NLsolve]] +deps = ["Distances", "LineSearches", "LinearAlgebra", "NLSolversBase", "Printf", "Reexport"] +git-tree-sha1 = "019f12e9a1a7880459d0173c182e6a99365d7ac1" +uuid = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" +version = "4.5.1" [[deps.NaNMath]] deps = ["OpenLibm_jll"] @@ -1634,9 +1657,9 @@ version = "1.2.0" [[deps.OrdinaryDiffEqBDF]] deps = ["ADTypes", "ArrayInterface", "DiffEqBase", "FastBroadcast", "LinearAlgebra", "MacroTools", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "OrdinaryDiffEqNonlinearSolve", "OrdinaryDiffEqSDIRK", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "StaticArrays", "TruncatedStacktraces"] -git-tree-sha1 = "42755bd13fe56e9d9ce1bc005f8b206a6b56b731" +git-tree-sha1 = "9124a686af119063bb4d3a8f87044a8f312fcad9" uuid = "6ad6398a-0878-4a85-9266-38940aa047c8" -version = "1.5.1" +version = "1.6.0" [[deps.OrdinaryDiffEqCore]] deps = ["ADTypes", "Accessors", "Adapt", "ArrayInterface", "DataStructures", "DiffEqBase", "DocStringExtensions", "EnumX", "FastBroadcast", "FastClosures", "FastPower", "FillArrays", "FunctionWrappersWrappers", "InteractiveUtils", "LinearAlgebra", "Logging", "MacroTools", "MuladdMacro", "Polyester", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "SciMLStructures", "SimpleUnPack", "Static", "StaticArrayInterface", "StaticArraysCore", "SymbolicIndexingInterface", "TruncatedStacktraces"] @@ -1656,9 +1679,9 @@ version = "1.4.0" [[deps.OrdinaryDiffEqDifferentiation]] deps = ["ADTypes", "ArrayInterface", "ConcreteStructs", "ConstructionBase", "DiffEqBase", "DifferentiationInterface", "FastBroadcast", "FiniteDiff", "ForwardDiff", "FunctionWrappersWrappers", "LinearAlgebra", "LinearSolve", "OrdinaryDiffEqCore", "SciMLBase", "SciMLOperators", "SparseArrays", "SparseMatrixColorings", "StaticArrayInterface", "StaticArrays"] -git-tree-sha1 = "c78060115fa4ea9d70ac47fa49496acbc630aefa" +git-tree-sha1 = "efecf0c4cc44e16251b0e718f08b0876b2a82b80" uuid = "4302a76b-040a-498a-8c04-15b101fed76b" -version = "1.9.1" +version = "1.10.0" [[deps.OrdinaryDiffEqExplicitRK]] deps = ["DiffEqBase", "FastBroadcast", "LinearAlgebra", "MuladdMacro", "OrdinaryDiffEqCore", "RecursiveArrayTools", "Reexport", "TruncatedStacktraces"] @@ -1764,9 +1787,9 @@ version = "1.1.0" [[deps.OrdinaryDiffEqRosenbrock]] deps = ["ADTypes", "DiffEqBase", "DifferentiationInterface", "FastBroadcast", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "LinearSolve", "MacroTools", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "Polyester", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "Static"] -git-tree-sha1 = "063e5ff1447b3869856ed264b6dcbb21e6e8bdb0" +git-tree-sha1 = "1ce0096d920e95773220e818f29bf4b37ea2bb78" uuid = "43230ef6-c299-4910-a778-202eb28ce4ce" -version = "1.10.1" +version = "1.11.0" [[deps.OrdinaryDiffEqSDIRK]] deps = ["ADTypes", "DiffEqBase", "FastBroadcast", "LinearAlgebra", "MacroTools", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "OrdinaryDiffEqNonlinearSolve", "RecursiveArrayTools", "Reexport", "SciMLBase", "TruncatedStacktraces"] @@ -2128,9 +2151,9 @@ version = "0.1.0" [[deps.SciMLBase]] deps = ["ADTypes", "Accessors", "Adapt", "ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "EnumX", "FunctionWrappersWrappers", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "Markdown", "Moshi", "PrecompileTools", "Preferences", "Printf", "RecipesBase", "RecursiveArrayTools", "Reexport", "RuntimeGeneratedFunctions", "SciMLOperators", "SciMLStructures", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface"] -git-tree-sha1 = "ce947672206f6a3a2bee1017c690cfd5fd82d897" +git-tree-sha1 = "d7bef263e23c7f5392071e4538b1949cee7ae458" uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462" -version = "2.96.0" +version = "2.99.0" [deps.SciMLBase.extensions] SciMLBaseChainRulesCoreExt = "ChainRulesCore" @@ -2357,6 +2380,12 @@ weakdeps = ["ChainRulesCore", "InverseFunctions"] StatsFunsChainRulesCoreExt = "ChainRulesCore" StatsFunsInverseFunctionsExt = "InverseFunctions" +[[deps.StochasticDiffEq]] +deps = ["ADTypes", "Adapt", "ArrayInterface", "DataStructures", "DiffEqBase", "DiffEqNoiseProcess", "DocStringExtensions", "FastPower", "FiniteDiff", "ForwardDiff", "JumpProcesses", "LevyArea", "LinearAlgebra", "Logging", "MuladdMacro", "NLsolve", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "OrdinaryDiffEqNonlinearSolve", "Random", "RandomNumbers", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "SparseArrays", "StaticArrays", "UnPack"] +git-tree-sha1 = "2992af2739fdcf5862b12dcf53a5f6e3e4acd358" +uuid = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0" +version = "6.80.0" + [[deps.StrideArraysCore]] deps = ["ArrayInterface", "CloseOpenIntervals", "IfElse", "LayoutPointers", "LinearAlgebra", "ManualMemory", "SIMDTypes", "Static", "StaticArrayInterface", "ThreadingUtilities"] git-tree-sha1 = "f35f6ab602df8413a50c4a25ca14de821e8605fb" @@ -2463,9 +2492,9 @@ version = "1.0.1" [[deps.Tables]] deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "OrderedCollections", "TableTraits"] -git-tree-sha1 = "598cd7c1f68d1e205689b1c2fe65a9f85846f297" +git-tree-sha1 = "f2c1efbc8f3a609aadf318094f8fc5204bdaf344" uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" -version = "1.12.0" +version = "1.12.1" [[deps.Tar]] deps = ["ArgTools", "SHA"] @@ -2556,14 +2585,16 @@ version = "0.4.1" [[deps.Unitful]] deps = ["Dates", "LinearAlgebra", "Random"] -git-tree-sha1 = "d62610ec45e4efeabf7032d67de2ffdea8344bed" +git-tree-sha1 = "02c1ac8104c9cf941395db79c611483909c04c7d" uuid = "1986cc42-f94f-5a68-af5c-568840ba703d" -version = "1.22.1" -weakdeps = ["ConstructionBase", "InverseFunctions"] +version = "1.23.0" +weakdeps = ["ConstructionBase", "ForwardDiff", "InverseFunctions", "Printf"] [deps.Unitful.extensions] ConstructionBaseUnitfulExt = "ConstructionBase" + ForwardDiffExt = "ForwardDiff" InverseFunctionsUnitfulExt = "InverseFunctions" + PrintfExt = "Printf" [[deps.UnitfulLatexify]] deps = ["LaTeXStrings", "Latexify", "Unitful"] @@ -2725,28 +2756,28 @@ uuid = "12413925-8142-5f55-bb0e-6d7ca50bb09b" version = "0.4.0+1" [[deps.Xorg_xcb_util_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll"] -git-tree-sha1 = "e7fd7b2881fa2eaa72717420894d3938177862d1" +deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libxcb_jll"] +git-tree-sha1 = "68da27247e7d8d8dafd1fcf0c3654ad6506f5f97" uuid = "2def613f-5ad1-5310-b15b-b15d46f528f5" -version = "0.4.0+1" +version = "0.4.1+0" [[deps.Xorg_xcb_util_keysyms_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"] -git-tree-sha1 = "d1151e2c45a544f32441a567d1690e701ec89b00" +deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_xcb_util_jll"] +git-tree-sha1 = "44ec54b0e2acd408b0fb361e1e9244c60c9c3dd4" uuid = "975044d2-76e6-5fbe-bf08-97ce7c6574c7" -version = "0.4.0+1" +version = "0.4.1+0" [[deps.Xorg_xcb_util_renderutil_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"] -git-tree-sha1 = "dfd7a8f38d4613b6a575253b3174dd991ca6183e" +deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_xcb_util_jll"] +git-tree-sha1 = "5b0263b6d080716a02544c55fdff2c8d7f9a16a0" uuid = "0d47668e-0667-5a69-a72c-f761630bfb7e" -version = "0.3.9+1" +version = "0.3.10+0" [[deps.Xorg_xcb_util_wm_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"] -git-tree-sha1 = "e78d10aab01a4a154142c5006ed44fd9e8e31b67" +deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_xcb_util_jll"] +git-tree-sha1 = "f233c83cad1fa0e70b7771e0e21b061a116f2763" uuid = "c22f9ab0-d5fe-5066-847c-f4bb1cd4e361" -version = "0.4.1+1" +version = "0.4.2+0" [[deps.Xorg_xkbcomp_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libxkbfile_jll"] @@ -2837,10 +2868,10 @@ uuid = "1183f4f0-6f2a-5f1a-908b-139f9cdfea6f" version = "0.2.2+0" [[deps.libevdev_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "141fe65dc3efabb0b1d5ba74e91f6ad26f84cc22" +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "56d643b57b188d30cccc25e331d416d3d358e557" uuid = "2db6ffa8-e38f-5e21-84af-90c45d0032cc" -version = "1.11.0+0" +version = "1.13.4+0" [[deps.libfdk_aac_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] @@ -2873,10 +2904,10 @@ uuid = "f27f6e37-5d2b-51aa-960f-b287f2bc3b7a" version = "1.3.7+2" [[deps.mtdev_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "814e154bdb7be91d78b6802843f76b6ece642f11" +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "b4d631fd51f2e9cdd93724ae25b2efc198b059b1" uuid = "009596ad-96f7-51b1-9f1b-5ce2d5e8a71e" -version = "1.1.6+0" +version = "1.1.7+0" [[deps.nghttp2_jll]] deps = ["Artifacts", "Libdl"] diff --git a/benchmarks/Jumps/MultivariateHawkes.jmd b/benchmarks/HybridJumps/MultivariateHawkes.jmd similarity index 88% rename from benchmarks/Jumps/MultivariateHawkes.jmd rename to benchmarks/HybridJumps/MultivariateHawkes.jmd index ecf52e137..4fd8be9c8 100644 --- a/benchmarks/Jumps/MultivariateHawkes.jmd +++ b/benchmarks/HybridJumps/MultivariateHawkes.jmd @@ -117,6 +117,7 @@ end function hawkes_problem( p, agg; + vr_agg = VR_FRM(), u = [0.0], tspan = (0.0, 50.0), save_positions = (false, true), @@ -125,7 +126,7 @@ function hawkes_problem( ) oprob = ODEProblem(f!, u, tspan, p) jumps = hawkes_jump(u, g; use_recursion) - jprob = JumpProblem(oprob, agg, jumps...; save_positions = save_positions) + jprob = JumpProblem(oprob, agg, jumps...; vr_aggregator = vr_agg, save_positions = save_positions) return jprob end ``` @@ -840,6 +841,100 @@ let fig = plot( end ``` +# Benchmarking Variable Rate Aggregators + +We benchmark the variable rate aggregators (`VR_Direct`, `VR_DirectFW`, `VR_FRM`) for the Multivariate Hawkes process, using the same setup as above: networks from `1` to `50` nodes, `tspan=(0.0, 25.0)`, `\lambda=0.5`, `\alpha=0.1`, `\beta=5.0`, and 50 trajectories with a 10-second limit per configuration. We test both recursive and brute-force formulations. + +```julia +vr_aggs = [ + (VR_Direct(), Tsit5(), false, "VR_Direct (brute-force)"), + (VR_DirectFW(), Tsit5(), false, "VR_DirectFW (brute-force)"), + (VR_FRM(), Tsit5(), false, "VR_FRM (brute-force)"), + (VR_Direct(), Tsit5(), true, "VR_Direct (recursive)"), + (VR_DirectFW(), Tsit5(), true, "VR_DirectFW (recursive)"), + (VR_FRM(), Tsit5(), true, "VR_FRM (recursive)"), +] + +tspan = (0.0, 25.0) +p = (0.5, 0.1, 5.0) +Vs = append!([1], 5:5:95) +Gs = [erdos_renyi(V, 0.2, seed = 6221) for V in Vs] + +vr_bs = Vector{Vector{BenchmarkTools.Trial}}() + +for (vr_agg, stepper, use_recursion, label) in vr_aggs + @info label + global _stepper = stepper + push!(vr_bs, Vector{BenchmarkTools.Trial}()) + _vr_bs = vr_bs[end] + for (i, G) in enumerate(Gs) + local g = [neighbors(G, i) for i in 1:nv(G)] + local u = [0.0 for i in 1:nv(G)] + if use_recursion + global h = zeros(eltype(tspan), nv(G)) + global urate = zeros(eltype(u), nv(G)) + global ϕ = zeros(eltype(tspan), nv(G)) + _p = (p[1], p[2], p[3], h, urate, ϕ) + else + global h = [eltype(tspan)[] for _ in 1:nv(G)] + global urate = zeros(eltype(u), nv(G)) + _p = (p[1], p[2], p[3], h, urate) + end + global jump_prob = hawkes_problem(_p, Direct(); vr_agg, u, tspan, g, use_recursion) + trial = try + if use_recursion + @benchmark( + solve(jump_prob, _stepper), + setup = (h .= 0; urate .= 0; ϕ .= 0), + samples = 50, + evals = 1, + seconds = 10, + ) + else + @benchmark( + solve(jump_prob, _stepper), + setup = ([empty!(_h) for _h in h]; urate .= 0), + samples = 50, + evals = 1, + seconds = 10, + ) + end + catch e + BenchmarkTools.Trial( + BenchmarkTools.Parameters(samples=50, evals=1, seconds=10), + ) + end + push!(_vr_bs, trial) + if (nv(G) == 1 || nv(G) % 10 == 0) + median_time = + length(trial) > 0 ? "$(BenchmarkTools.prettytime(median(trial.times)))" : "nan" + println("algo=$label, V=$(nv(G)), length=$(length(trial.times)), median time=$median_time") + end + end +end +``` + +```julia +let fig = plot( + yscale = :log10, + xlabel = "V", + ylabel = "Time (ns)", + legend_position = :outertopright, +) + for (i, (vr_agg, _, use_recursion, label)) in enumerate(vr_aggs) + _bs, _Vs = [], [] + for (j, b) in enumerate(vr_bs[i]) + if length(b) == 50 + push!(_bs, median(b.times)) + push!(_Vs, Vs[j]) + end + end + plot!(_Vs, _bs, label=label) + end + title!("Variable Rate Simulations, 50 samples: nodes × time") +end +``` + # References [1] D. J. Daley and D. Vere-Jones. An Introduction to the Theory of Point Processes: Volume I: Elementary Theory and Methods. Probability and Its Applications, An Introduction to the Theory of Point Processes. Springer-Verlag, 2 edition. doi:10.1007/b97277. diff --git a/benchmarks/HybridJumps/Project.toml b/benchmarks/HybridJumps/Project.toml index 6838c8426..f3a48d133 100644 --- a/benchmarks/HybridJumps/Project.toml +++ b/benchmarks/HybridJumps/Project.toml @@ -2,11 +2,14 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" +DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" JumpProcesses = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5" LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" PiecewiseDeterministicMarkovProcesses = "86206cdf-4603-54e0-bd58-22a2dcbf57aa" @@ -14,7 +17,9 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" SciMLBenchmarks = "31c91b34-3c75-11e9-0341-95557aab0344" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" +StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0" [compat] BenchmarkTools = "1.4" diff --git a/benchmarks/HybridJumps/Synapse.jmd b/benchmarks/HybridJumps/Synapse.jmd index 66e21e280..4eee4ed9b 100644 --- a/benchmarks/HybridJumps/Synapse.jmd +++ b/benchmarks/HybridJumps/Synapse.jmd @@ -1640,7 +1640,7 @@ function J_synapse(p_synapse::SynapseParams, nu) # we order the jumps in their order they appear in the dependency graph jumps = JumpSet(; - constant_jumps = [ + constant_jumps = ( # AMPA #2line-GO @j_jump(1, p_synapse, nu, 4 * AMPA_k1 * p.Glu * p.xd[1]), # 1 @@ -1740,8 +1740,8 @@ function J_synapse(p_synapse::SynapseParams, nu) @j_jump(97, p_synapse, nu, GABA_r_c1 * p.xd[49]), # 97 @j_jump(98, p_synapse, nu, GABA_r_ro2 * p.xd[48]), # 98 @j_jump(99, p_synapse, nu, GABA_r_c2 * p.xd[50]), # 99 - ], - variable_jumps = [ + ), + variable_jumps = ( # R-type VGCC @j_jump( 56, @@ -1913,7 +1913,7 @@ function J_synapse(p_synapse::SynapseParams, nu) @j_jump(77, p_synapse, nu, p.xd[37] * P_rate, 1, typemax(Float64)), # 77 @j_jump(78, p_synapse, nu, p.xd[36] * P_rate, 1, typemax(Float64)), # 78 @j_jump(79, p_synapse, nu, p.xd[38] * D_rate, 1, typemax(Float64)), # 79 - ], + ), ) return jumps @@ -1955,6 +1955,7 @@ function SynapseProblem( nu, algo::T, agg = nothing; + vr_agg = VR_FRM(), saveat = [], save_everystep = isempty(saveat), kwargs..., @@ -2128,6 +2129,7 @@ function SynapseProblem( nu, algo, agg; + vr_agg = VR_FRM(), jumps = nothing, save_positions = (false, true), saveat = [], @@ -2158,6 +2160,7 @@ function SynapseProblem( oprob, agg, jumps; + vr_aggregator = vr_agg, dep_graph, save_positions, saveat, @@ -2186,6 +2189,7 @@ function evolveSynapse( nu, algos, agg = nothing; + vr_agg, progress = false, abstol = 1e-8, reltol = 1e-7, @@ -2205,6 +2209,7 @@ function evolveSynapse( nu, algos, agg; + vr_agg, progress, abstol, reltol, @@ -2228,6 +2233,7 @@ function evolveSynapse_noformat( nu, algos, agg = nothing; + vr_agg, progress = false, abstol = 1e-8, reltol = 1e-7, @@ -2274,6 +2280,7 @@ function evolveSynapse_noformat( nu, algos[1], agg; + vr_agg, jumps, reltol, abstol, @@ -2296,6 +2303,7 @@ function evolveSynapse_noformat( nu, algos[2], agg; + vr_agg, jumps, reltol, abstol, @@ -2415,12 +2423,35 @@ const algorithms = ( ( label = "PDMP", agg = nothing, + vr_agg = VR_FRM(), solver = (CHV(solver), CHV(solver)), saveat = [], ), ( label = "Coevolve", agg = Coevolve(), + vr_agg = VR_FRM(), + solver = (solver, solver), + saveat = 1 / p_synapse.sampling_rate, + ), + ( + label = "Direct(VR_FRM)", + agg = Direct(), + vr_agg = VR_FRM(), + solver = (solver, solver), + saveat = 1 / p_synapse.sampling_rate, + ), + ( + label = "Direct(VR_Direct)", + agg = Direct(), + vr_agg = VR_Direct(), + solver = (solver, solver), + saveat = 1 / p_synapse.sampling_rate, + ), + ( + label = "Direct(VR_DirectFW)", + agg = Direct(), + vr_agg = VR_DirectFW(), solver = (solver, solver), saveat = 1 / p_synapse.sampling_rate, ), @@ -2446,6 +2477,7 @@ for algo in algorithms nu, algo.solver, algo.agg; + algo.vr_agg, save_positions = (false, true), saveat = algo.saveat, save_everystep = false, @@ -2492,6 +2524,7 @@ for algo in algorithms nu, $(algo).solver, $(algo).agg; + $(algo).vr_agg, save_positions = (false, true), saveat = $(algo).saveat, save_everystep = false, diff --git a/benchmarks/HybridJumps/VR_Aggregator_Benchmark.jmd b/benchmarks/HybridJumps/VR_Aggregator_Benchmark.jmd new file mode 100644 index 000000000..a93e2fea5 --- /dev/null +++ b/benchmarks/HybridJumps/VR_Aggregator_Benchmark.jmd @@ -0,0 +1,477 @@ +--- +title: Benchmarking Variable Rate Aggregator +author: Siva Sathyaseelan D N, Chris Rackauckas, Samuel Isaacson +weave_options: + fig_ext : ".png" +--- + +```julia +using DiffEqBase, Catalyst, JumpProcesses, OrdinaryDiffEq, StochasticDiffEq +using Random, LinearSolve, StableRNGs, BenchmarkTools, Plots, LinearAlgebra +fmt = :png +width_px, height_px = default(:size) +rng = StableRNG(12345) +``` + +# Introduction + +This document benchmarks the performance of variable rate jumps in `JumpProcesses.jl` and visualizes example solution trajectories. The benchmark compares `VR_Direct`, `VR_DirectFW`, and `VR_FRM` aggregators for variable rate jumps, and includes a constant rate jump model with `Direct`. Visualization shows state variables vs. time to verify behavior. + +The test cases are: +1. **Scalar ODE with Variable Rate Jumps**: Solved with `Tsit5` and `Rosenbrock23` (autodiff). +2. **Complex ODE with Variable Rate Jump**: Solved with `Tsit5`. +3. **DNA Gene Model**: ODE with 10 variable rate jumps from the RSSA paper, solved with `Tsit5`. +4. **Negative Feedback Gene Expression**: Constant rate jumps solved with `SSAStepper`, variable rate jumps with `Tsit5`, from Marchetti et al. (2017). + +For visualization, we solve one trajectory per test case with 2 jumps (10 for Test 3, 8 for Test 4). For benchmarking, we vary jumps from 1 to 20 for Tests 1 and 2, use fixed jumps for Test 3 (10) and Test 4 (8), running 50 trajectories. Benchmarking saves only at the final time with `save_positions=(false, false)`. + +# Benchmark and Plot Test 1 + +We benchmark Test 1 for 1 to 20 jumps, running 50 trajectories, and plot mean execution times as a line plot. + +```julia +let + algorithms = Tuple{Any, Any, String, String}[ + (VR_Direct(), Tsit5(), "VR_Direct", "Test 1 Tsit5 (VR_Direct)"), + (VR_DirectFW(), Tsit5(), "VR_DirectFW", "Test 1 Tsit5 (VR_DirectFW)"), + (VR_FRM(), Tsit5(), "VR_FRM", "Test 1 Tsit5 (VR_FRM)"), + (VR_Direct(), Rosenbrock23(), "VR_Direct", "Test 1 Rosenbrock23 (autodiff, VR_Direct)"), + (VR_DirectFW(), Rosenbrock23(), "VR_DirectFW", "Test 1 Rosenbrock23 (autodiff, VR_DirectFW)"), + (VR_FRM(), Rosenbrock23(), "VR_FRM", "Test 1 Rosenbrock23 (autodiff, VR_FRM)"), + ] + + function create_test1_problem(num_jumps, vr_aggregator, solver) + f = (du, u, p, t) -> (du[1] = u[1]) + prob = ODEProblem(f, [0.2], (0.0, 10.0)) + jumps = [VariableRateJump((u, p, t) -> u[1], (integrator) -> (integrator.u[1] = integrator.u[1] / 2.0); interp_points=20) for _ in 1:num_jumps] + jump_prob = JumpProblem(prob, Direct(), jumps...; vr_aggregator=vr_aggregator, rng=rng, save_positions=(false, false)) + ensemble_prob = EnsembleProblem(jump_prob) + return ensemble_prob, jump_prob + end + + num_jumps_range = append!([1], 5:5:20) + bs = Vector{Vector{BenchmarkTools.Trial}}() + errors = Dict{String, Vector{String}}() + + for (algo, stepper, agg_name, label) in algorithms + @info "Benchmarking $label" + push!(bs, Vector{BenchmarkTools.Trial}()) + errors[label] = String[] + _bs = bs[end] + for var in num_jumps_range + ensemble_prob, jump_prob = create_test1_problem(var, algo, stepper) + trial = try + @benchmark( + solve($jump_prob, $stepper, saveat=[$jump_prob.prob.tspan[2]]), + samples=50, + evals=1, + seconds=100 + ) + catch e + push!(errors[label], "Error at Num Jumps = $var: $(sprint(showerror, e))") + BenchmarkTools.Trial(BenchmarkTools.Parameters(samples=50, evals=1, seconds=100)) + end + push!(_bs, trial) + mean_time = length(trial) > 0 ? "$(BenchmarkTools.prettytime(mean(trial.times)))" : "nan" + println("algo=$label, Num Jumps = $var, length = $(length(trial.times)), mean time = $mean_time") + end + end + + # Log errors + for (label, err_list) in errors + if !isempty(err_list) + @warn "Errors for $label:" + for err in err_list + println(err) + end + end + end + + # Plot results + fig = plot( + yscale=:log10, + xlabel="Number of Jumps", + ylabel="Time (ns)", + legend_position=:outertopright, + title="Test 1: Simulations, 50 samples" + ) + for (i, (algo, stepper, agg_name, label)) in enumerate(algorithms) + _bs, _vars = [], [] + for (j, b) in enumerate(bs[i]) + if length(b) == 50 + push!(_bs, mean(b.times)) + push!(_vars, num_jumps_range[j]) + end + end + if !isempty(_bs) + plot!(_vars, _bs, label=label) + else + @warn "No valid data for $label in Test 1" + end + end + display(plot(fig, layout=(1, 1), format=fmt, size=(width_px, height_px))) +end +``` + +# Benchmark and Plot Test 2 + +We benchmark Test 2 for 1 to 20 jumps, running 50 trajectories, and plot mean execution times as a line plot. + +```julia +let + algorithms = Tuple{Any, Any, String, String}[ + (VR_Direct(), Tsit5(), "VR_Direct", "Test 2 Tsit5 (VR_Direct)"), + (VR_DirectFW(), Tsit5(), "VR_DirectFW", "Test 2 Tsit5 (VR_DirectFW)"), + (VR_FRM(), Tsit5(), "VR_FRM", "Test 2 Tsit5 (VR_FRM)"), + ] + + function create_test2_problem(num_jumps, vr_aggregator, solver) + f4 = (dx, x, p, t) -> (dx[1] = x[1]) + rate4 = (x, p, t) -> t + affect4! = (integrator) -> (integrator.u[1] = integrator.u[1] * 0.5) + prob = ODEProblem(f4, [1.0 + 0.0im], (0.0, 6.0)) + jumps = [VariableRateJump(rate4, affect4!) for _ in 1:num_jumps] + jump_prob = JumpProblem(prob, Direct(), jumps...; vr_aggregator=vr_aggregator, rng=rng, save_positions=(false, false)) + ensemble_prob = EnsembleProblem(jump_prob) + return ensemble_prob, jump_prob + end + + num_jumps_range = append!([1], 5:5:20) + bs = Vector{Vector{BenchmarkTools.Trial}}() + errors = Dict{String, Vector{String}}() + + for (algo, stepper, agg_name, label) in algorithms + @info "Benchmarking $label" + push!(bs, Vector{BenchmarkTools.Trial}()) + errors[label] = String[] + _bs = bs[end] + for var in num_jumps_range + ensemble_prob, jump_prob = create_test2_problem(var, algo, stepper) + trial = try + @benchmark( + solve($jump_prob, $stepper, saveat=[$jump_prob.prob.tspan[2]]), + samples=50, + evals=1, + seconds=100 + ) + catch e + push!(errors[label], "Error at Num Jumps = $var: $(sprint(showerror, e))") + BenchmarkTools.Trial(BenchmarkTools.Parameters(samples=50, evals=1, seconds=100)) + end + push!(_bs, trial) + mean_time = length(trial) > 0 ? "$(BenchmarkTools.prettytime(mean(trial.times)))" : "nan" + println("algo=$label, Num Jumps = $var, length = $(length(trial.times)), mean time = $mean_time") + end + end + + # Log errors + for (label, err_list) in errors + if !isempty(err_list) + @warn "Errors for $label:" + for err in err_list + println(err) + end + end + end + + # Plot results + fig = plot( + yscale=:log10, + xlabel="Number of Jumps", + ylabel="Time (ns)", + legend_position=:outertopright, + title="Test 2: Simulations, 50 samples" + ) + for (i, (algo, stepper, agg_name, label)) in enumerate(algorithms) + _bs, _vars = [], [] + for (j, b) in enumerate(bs[i]) + if length(b) == 50 + push!(_bs, mean(b.times)) + push!(_vars, num_jumps_range[j]) + end + end + if !isempty(_bs) + plot!(_vars, _bs, label=label) + else + @warn "No valid data for $label in Test 2" + end + end + display(plot(fig, layout=(1, 1), format=fmt, size=(width_px, height_px))) +end +``` + +# Benchmark and Plot Test 3 + +We benchmark Test 3 for fixed 10 jumps, running 50 trajectories, and plot mean execution times as a bar plot. + +```julia +let + algorithms = Tuple{Any, Any, String, String}[ + (VR_Direct(), Tsit5(), "VR_Direct", "Test 3 Tsit5 (VR_Direct, DNA Model)"), + (VR_DirectFW(), Tsit5(), "VR_DirectFW", "Test 3 Tsit5 (VR_DirectFW, DNA Model)"), + (VR_FRM(), Tsit5(), "VR_FRM", "Test 3 Tsit5 (VR_FRM, DNA Model)"), + ] + + function create_test3_problem(num_jumps, vr_aggregator, solver) + r = [0.043, 0.0007, 0.0715, 0.0039, 0.0199, 0.4791, 0.00019, 0.8765, 0.083, 0.5] + k = -log(2) / 30 + u0 = [10.0, 10.0, 30.0, 0.0, 0.0, 0.0] # [DNA, M, D, RNA, DNAD, DNA2D] + tspan = (0.0, 120.0) + + function f_dna(du, u, p, t) + du .= 0.0 + nothing + end + + function rate1(u, p, t) r[1] * u[4] end + function affect1!(integrator) integrator.u[2] += 1; nothing end + jump1 = VariableRateJump(rate1, affect1!) + + function rate2(u, p, t) r[2] * u[2] end + function affect2!(integrator) integrator.u[2] -= 1; nothing end + jump2 = VariableRateJump(rate2, affect2!) + + function rate3(u, p, t) r[3] * u[5] end + function affect3!(integrator) integrator.u[4] += 1; nothing end + jump3 = VariableRateJump(rate3, affect3!) + + function rate4(u, p, t) r[4] * u[4] end + function affect4!(integrator) integrator.u[4] -= 1; nothing end + jump4 = VariableRateJump(rate4, affect4!) + + function rate5(u, p, t) r[5] * exp(k * t) * u[1] * u[3] end + function affect5!(integrator) integrator.u[1] -= 1; integrator.u[3] -= 1; integrator.u[5] += 1; nothing end + jump5 = VariableRateJump(rate5, affect5!) + + function rate6(u, p, t) r[6] * u[5] end + function affect6!(integrator) integrator.u[5] -= 1; integrator.u[1] += 1; integrator.u[3] += 1; nothing end + jump6 = VariableRateJump(rate6, affect6!) + + function rate7(u, p, t) r[7] * exp(k * t) * u[5] * u[3] end + function affect7!(integrator) integrator.u[5] -= 1; integrator.u[3] -= 1; integrator.u[6] += 1; nothing end + jump7 = VariableRateJump(rate7, affect7!) + + function rate8(u, p, t) r[8] * u[6] end + function affect8!(integrator) integrator.u[6] -= 1; integrator.u[1] += 1; integrator.u[3] += 1; nothing end + jump8 = VariableRateJump(rate8, affect8!) + + function rate9(u, p, t) r[9] * exp(k * t) * u[2] * (u[2] - 1) / 2 end + function affect9!(integrator) integrator.u[2] -= 2; integrator.u[3] += 1; nothing end + jump9 = VariableRateJump(rate9, affect9!) + + function rate10(u, p, t) r[10] * u[3] end + function affect10!(integrator) integrator.u[3] -= 1; integrator.u[2] += 2; nothing end + jump10 = VariableRateJump(rate10, affect10!) + + prob = ODEProblem(f_dna, u0, tspan) + jumps = (jump1, jump2, jump3, jump4, jump5, jump6, jump7, jump8, jump9, jump10) + jump_prob = JumpProblem(prob, Direct(), jumps...; vr_aggregator=vr_aggregator, rng=rng, save_positions=(false, false)) + ensemble_prob = EnsembleProblem(jump_prob) + return ensemble_prob, jump_prob + end + + num_jumps_range = [10] + bs = Vector{Vector{BenchmarkTools.Trial}}() + errors = Dict{String, Vector{String}}() + + for (algo, stepper, agg_name, label) in algorithms + @info "Benchmarking $label" + push!(bs, Vector{BenchmarkTools.Trial}()) + errors[label] = String[] + _bs = bs[end] + for var in num_jumps_range + ensemble_prob, jump_prob = create_test3_problem(var, algo, stepper) + trial = try + @benchmark( + solve($jump_prob, $stepper, saveat=[$jump_prob.prob.tspan[2]]), + samples=50, + evals=1, + seconds=100 + ) + catch e + push!(errors[label], "Error at Num Jumps = $var: $(sprint(showerror, e))") + BenchmarkTools.Trial(BenchmarkTools.Parameters(samples=50, evals=1, seconds=100)) + end + push!(_bs, trial) + mean_time = length(trial) > 0 ? "$(BenchmarkTools.prettytime(mean(trial.times)))" : "nan" + println("algo=$label, Num Jumps = $var, length = $(length(trial.times)), mean time = $mean_time") + end + end + + # Log errors + for (label, err_list) in errors + if !isempty(err_list) + @warn "Errors for $label:" + for err in err_list + println(err) + end + end + end + + # Plot results + fig = plot( + yscale=:log10, + xlabel="Method", + ylabel="Time (ns)", + title="Test 3: DNA Model, 10 Jumps, 50 samples", + xticks=(1:length(algorithms), [split(a[3], " (")[1] for a in algorithms]), + xrotation=45 + ) + means = [] + for (i, (algo, stepper, agg_name, label)) in enumerate(algorithms) + b = bs[i][1] # Single jump count (10) + if length(b) == 50 + push!(means, mean(b.times)) + else + push!(means, NaN) + @warn "No valid data for $label in Test 3" + end + end + bar!(1:length(algorithms), means, label="", fillalpha=0.7) + display(plot(fig, layout=(1, 1), format=fmt, size=(width_px, height_px))) +end +``` + +# Benchmark and Plot Test 4 + +We benchmark Test 4 for fixed 8 jumps, running 50 trajectories, and plot mean execution times as a bar plot. + +```julia +let + algorithms = Tuple{Any, Any, String, String}[ + #(Direct(), SSAStepper(), "Direct", "Test 4 SSAStepper (Direct, NegFeedback, Constant Rate)"), + (VR_Direct(), Tsit5(), "VR_Direct", "Test 4 Tsit5 (VR_Direct, NegFeedback, Variable Rate)"), + (VR_DirectFW(), Tsit5(), "VR_DirectFW", "Test 4 Tsit5 (VR_DirectFW, NegFeedback, Variable Rate)"), + (VR_FRM(), Tsit5(), "VR_FRM", "Test 4 Tsit5 (VR_FRM, NegFeedback, Variable Rate)"), + ] + + function create_test4_problem(num_jumps, aggregator, solver) + rn = @reaction_network begin + c1, G --> G + M + c2, M --> M + P + c3, M --> 0 + c4, P --> 0 + c5, 2P --> P2 + c6, P2 --> 2P + c7, P2 + G --> P2G + c8, P2G --> P2 + G + end + rnpar = [:c1 => 0.09, :c2 => 0.05, :c3 => 0.001, :c4 => 0.0009, :c5 => 0.00001, + :c6 => 0.0005, :c7 => 0.005, :c8 => 0.9] + u0 = [:G => 500, :M => 0, :P => 0, :P2 => 0, :P2G => 0] + tspan = (0.0, 100.0) + + if aggregator isa Direct + prob = DiscreteProblem(rn, u0, tspan, rnpar) + jump_prob = JumpProblem(rn, prob, Direct(), rng=rng, save_positions=(false, false)) + ensemble_prob = EnsembleProblem(jump_prob) + else + function f_gene(du, u, p, t) + du .= 0.0 + nothing + end + u0_numeric = [500.0, 0.0, 0.0, 0.0, 0.0] + p_numeric = [0.09, 0.05, 0.001, 0.0009, 0.00001, 0.0005, 0.005, 0.9] + var_rate = true + prob = ODEProblem(f_gene, u0_numeric, tspan, (p_numeric, var_rate)) + + function rate1(u, p, t) p[1][1] * (p[2] ? (1 + 0.1 * sin(t)) : 1.0) * u[1] end + function affect1!(integrator) integrator.u[2] += 1; nothing end + jump1 = VariableRateJump(rate1, affect1!) + + function rate2(u, p, t) p[1][2] * u[2] end + function affect2!(integrator) integrator.u[3] += 1; nothing end + jump2 = VariableRateJump(rate2, affect2!) + + function rate3(u, p, t) p[1][3] * u[2] end + function affect3!(integrator) integrator.u[2] -= 1; nothing end + jump3 = VariableRateJump(rate3, affect3!) + + function rate4(u, p, t) p[1][4] * u[3] end + function affect4!(integrator) integrator.u[3] -= 1; nothing end + jump4 = VariableRateJump(rate4, affect4!) + + function rate5(u, p, t) p[1][5] * (p[2] ? (1 + 0.1 * cos(t)) : 1.0) * u[3] * (u[3] - 1) / 2 end + function affect5!(integrator) integrator.u[3] -= 2; integrator.u[4] += 1; nothing end + jump5 = VariableRateJump(rate5, affect5!) + + function rate6(u, p, t) p[1][6] * u[4] end + function affect6!(integrator) integrator.u[4] -= 1; integrator.u[3] += 2; nothing end + jump6 = VariableRateJump(rate6, affect6!) + + function rate7(u, p, t) p[1][7] * u[4] * u[1] end + function affect7!(integrator) integrator.u[4] -= 1; integrator.u[1] -= 1; integrator.u[5] += 1; nothing end + jump7 = VariableRateJump(rate7, affect7!) + + function rate8(u, p, t) p[1][8] * u[5] end + function affect8!(integrator) integrator.u[5] -= 1; integrator.u[4] += 1; integrator.u[1] += 1; nothing end + jump8 = VariableRateJump(rate8, affect8!) + + jumps = (jump1, jump2, jump3, jump4, jump5, jump6, jump7, jump8) + jump_prob = JumpProblem(prob, Direct(), jumps...; vr_aggregator=aggregator, rng=rng, save_positions=(false, false)) + ensemble_prob = EnsembleProblem(jump_prob) + end + return ensemble_prob, jump_prob + end + + num_jumps_range = [8] + bs = Vector{Vector{BenchmarkTools.Trial}}() + errors = Dict{String, Vector{String}}() + + for (algo, stepper, agg_name, label) in algorithms + @info "Benchmarking $label" + push!(bs, Vector{BenchmarkTools.Trial}()) + errors[label] = String[] + _bs = bs[end] + for var in num_jumps_range + ensemble_prob, jump_prob = create_test4_problem(var, algo, stepper) + trial = try + @benchmark( + solve($jump_prob, $stepper, saveat=[$jump_prob.prob.tspan[2]]), + samples=50, + evals=1, + seconds=1000 + ) + catch e + push!(errors[label], "Error at Num Jumps = $var: $(sprint(showerror, e))") + BenchmarkTools.Trial(BenchmarkTools.Parameters(samples=10, evals=1, seconds=1000)) + end + push!(_bs, trial) + mean_time = length(trial) > 0 ? "$(BenchmarkTools.prettytime(mean(trial.times)))" : "nan" + println("algo=$label, Num Jumps = $var, length = $(length(trial.times)), mean time = $mean_time") + end + end + + # Log errors + for (label, err_list) in errors + if !isempty(err_list) + @warn "Errors for $label:" + for err in err_list + println(err) + end + end + end + + # Plot results + fig = plot( + yscale=:log10, + xlabel="Method", + ylabel="Time (ns)", + title="Test 4: NegFeedback, 8 Jumps, 50 samples", + xticks=(1:length(algorithms), [split(a[3], " (")[1] for a in algorithms]), + xrotation=45 + ) + means = [] + for (i, (algo, stepper, agg_name, label)) in enumerate(algorithms) + b = bs[i][1] # Single jump count (8) + if length(b) == 50 + push!(means, mean(b.times)) + else + push!(means, NaN) + @warn "No valid data for $label in Test 4" + end + end + bar!(1:length(algorithms), means, label="", fillalpha=0.7) + display(plot(fig, layout=(1, 1), format=fmt, size=(width_px, height_px))) +end +``` \ No newline at end of file