SciML / JumpProcesses.jl

Build and simulate jump equations like Gillespie simulations and jump diffusions with constant and state-dependent rates and mix with differential equations and scientific machine learning (SciML)
https://docs.sciml.ai/JumpProcesses/stable/
Other
140 stars 35 forks source link

Choose an SSA if no SSA is passed in `JumpProblem`. #351

Closed Vilin97 closed 1 month ago

Vilin97 commented 1 year ago

If JumpProblem(prob, jumps) is called, an SSA is selected via the function below. The assumption is that Direct is fastest for a small network, while RSSACR is fastest for other networks. If RSSACR cannot be used because no species-to-jumps graph is given, then DirectCR is used instead.

# 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)

    # detect if a spatial SSA should be used
    !isnothing(spatial_system) && !isnothing(hopping_constants) && return DirectCRDirect

    # if variable rate jumps are present, return the only SSA that supports them
    num_vrjs(jumps) != 0 && return Coevolve

    # if the number of jumps is small, return the Direct
    num_majumps(jumps) + num_crjs(jumps) < 10 && return Direct

    # if there are only massaction jumps, we can any build dependency graph, so return the fastest aggregator
    num_crjs(jumps) == 0 && num_vrjs(jumps) == 0 && return RSSACR

    # in the presence of constant rate jumps, we rely on the given dependency graphs
    if !isnothing(vartojumps_map) && !isnothing(jumptovars_map)
        return RSSACR
    elseif !isnothing(dep_graph)
        return DirectCR
    else
        return Direct
    end
end
github-actions[bot] commented 1 year ago

Pull Request Test Coverage Report for Build 7514454839

Warning: This coverage report may be inaccurate.

We've detected an issue with your CI configuration that might affect the accuracy of this pull request's coverage report. To ensure accuracy in future PRs, please see these guidelines. A quick fix for this PR: rebase it; your next report should be accurate.


Totals Coverage Status
Change from base Build 6269620671: -1.1%
Covered Lines: 2070
Relevant Lines: 2332

💛 - Coveralls
ChrisRackauckas commented 9 months ago

Perfect is always the enemy of good for this kind of thing. I can tell you that 90% of users just use Gillespie's method if we document Direct as Gillespie's method. Even the simplest defaulting scheme can be merged and do better than the vast majority of users, so the details can always evolve but I think the biggest thing is to just get something with the right broad strokes and also update all of the tutorials to not hardcode a method. The latter is BTW very important, because there's still people who use DiffEq who just throw Tsit5 everywhere instead of using the default, and you really want to stop that habit via tutorials otherwise people will overuse Gillespie.

isaacsas commented 1 month ago

@Vilin97 will merge when I get tests passing.

@TorkelE expanded this a bit to allow not needing to pass SSAStepper too, i.e. with this PR this should now work

jprob = JumpProblem(rn, dprob)
sol = solve(jprob)
isaacsas commented 1 month ago

@TorkelE actually MTK and Catalyst need appropriate dispatches for JumpProblem to make this work, so not quite there yet for being usable in Catalyst.

isaacsas commented 1 month ago

@ChrisRackauckas if you could give a quick look to ensure I haven't added a dispatch better handled elsewhere (or shadowing one elsewhere) that would be helpful.

isaacsas commented 1 month ago

@Vilin97 thanks for your work on this and sorry for the long long delay!

Vilin97 commented 1 month ago

Thank you for cleaning up and adding a default solver too. Now we need to update the docs, right?

isaacsas commented 1 month ago

Yes, the docs need updates in places now. Using the defaults in the most introductory places, but explaining how to select solvers and such as they currently do in more advanced places.