-
-
Notifications
You must be signed in to change notification settings - Fork 34
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Large SSA support #16
Comments
For SSAs that require more information, maybe another option would be to just have the Also, I don't think large SSAs necessarily need to reorder the vector of jump rates or affects; one could always internally maintain a permutation vector (ignoring if this has performance implications). However, they do need to be able to only update a small subset of jump rates and species populations each step. You want to avoid loops over all reactions or all species during each step. |
With regard to handling larger networks in the current codebase; I’ve now played around with using FunctionWrappers for both Direct and FRM with a larger network. (1D diffusion through hopping with equal rates between 50 lattice sites, with 10 molecules per site initially -- this gives 98 total reactions). All I did was replace the tuple for It definitely makes a big difference over the tuple approach, but does seem to be a bit slower than BioSimulator. Benchmarks: (key Direct2 = regular Direct with tuples, DirectVec = Direct with function wrappers, FRM = FRM with tuples, FRMVec = FRM with function wrappers) Full solution paths saving: Solving with method: DiffEqJumpExtensions.Direct2, using SSAStepper Saving at fixed times: Solving with method: DiffEqJumpExtensions.Direct2, using SSAStepper and 25000 points I'd also like to try an approach like BioSimulator and Bionetgen, where there is just a global rate function that takes the network info and dynamically calculates a requested rate. I don't know why this would be faster for a smaller network, but perhaps for a system with thousands (or tens of thousands) of reactions it could be quicker to evaluate a global function with different parameters (the rate id) than to lookup and evaluate a dedicated function for each specific reaction... |
What did you leave
Oh you mean like a Did you profile to see where all of the time is spent? I noticed for instance that this is missing an |
I wasn't quite sure what type signature to use in
setting them up in the JumpProblem constructor. (Though perhaps I should replace the Float64 with the type of I haven't tried using the profiler yet, but that is on the TODO list... |
It should just ignore the output, so it would be helpful to do wrapped_f = x->(affect!(x);nothing) and then you can type the return as |
What type should the input have? Or should I just use
And yes, something like you write for stoich_rate. Though it might make sense to manually calculate the product, and the function would need to handle things like homodimerization reactions for which we use
|
You can use
Yup, that makes sense. |
Just to record my final benchmarking results: I'm not seeing a big difference in benchmarking the Direct method with using any of a vector of I put together a jump aggregator that uses a single global rate function and single global affect function for mass action reactions (similar to the BioSimulator design, but using vectors of vectors of pairs to store stochiometry instead dense/sparse matrices). For the direct method, on a 1D diffusion example with 510 reactions I find (DirectMA = mass action version, DirectVEC = function wrappers version) tf = 20 seconds, N = 256 lattice sites tf = 200 seconds, N = 256 lattice sites So the general conclusion seems to be that for generic jumps/affects, it would be good to have Direct.jl transition to FunctionWrappers once the system has enough jumps. For mass action systems it makes sense to have a specialized approach that uses a single function for rates and a single function for affects. |
How's the implementation? Is it easy to abstract and use the two jump types in the same SSA? That would mean the interface gains some complexity, but at least if the backend is independent of how the rates are formed then it's still not too bad. |
I don't think that would be tough. Right now I'm not using the You can see how I set it up for just mass action at: https://github.com/isaacsas/DiffEqJumpExtensions The relevant files are
For the stochiometry representation I tried a couple interfaces. A vector of vectors of pairs seemed to work well. (Here a given pair maps the id of a species to the stochiometric coefficient of the species for a given reaction.) The main issue I'm having is to figure out where the extra reaction info that the mass action routines need is passed in. An extra jump type would work, but what happens then for SSAs that want even more info (like a dependency graph)? Do we extend the jump types again, or for everything beyond stochiometry just pass this extra info through the aggregator algorithm or associated |
I'm not sure about the best way to do this. I think it may be impossible to plan that far ahead, and instead the code should be written for the current methods and keep expanding as needed. When we need it we'll probably think of the best solution, but trying to plan for it is really really hard since we don't even know what we need. |
Fair enough; right now I'm planning to just implement a mass action jump type (like the regular jump type), which will be another field in the |
We need to work out how we want to do large SSAs. We can build hybrid algorithms over the
RegularJump
s, but @isaacsas mentioned that large SSAs want to be able to do re-ordering so that's not the right data structure. But our current setup with tuples is optimized for less than 16 jumps, so if we want to useConstantRateJump
for this then we need to somehow useFunctionWrappers
and allow putting them in an array. However, there might be more information needed, in which case it might be useful to write SSAs on yet another jump type.The text was updated successfully, but these errors were encountered: