Example: Measuring activation time distribution #1210
Conversation
docs/Project.toml
Outdated
| @@ -40,6 +40,7 @@ SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" | |||
| SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1" | |||
| SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" | |||
| StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" | |||
| StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd" | |||
There was a problem hiding this comment.
A bit annoying to add a dependency for a tutorial. However, the density function feels quite important for stochastic -related analysis, so probably something users should know about and which we might want to use in other tutorials as well. Also, it is part of the Plots system, which should make it less of an addition.
There was a problem hiding this comment.
We probably should at some point just switch everything to Makie-verse given we are now pulling it in, but that is a longer-term goal I guess.
docs/src/model_simulation/examples/activation_time_distribution_measurement.md
Outdated
Show resolved
Hide resolved
| condition(u, t, integrator) = integrator[:X] > 100.0 | ||
| affect!(integrator) = terminate!(integrator) | ||
| callback = DiscreteCallback(condition, affect!) |
There was a problem hiding this comment.
Can you not use a functional affect and store this in the system?
There was a problem hiding this comment.
Does system events store non-trivial affects like termiante? If it is easier to do I am happy to change, however, it doesn't hurt to show how more complicated events can be implemented through callbacks.
| ```@example activation_time_distribution_measurement | ||
| function output_func(sol, i) | ||
| (sol.t[end] == tspan[2]) && @warn "A simulation did not activate during the given time span." | ||
| return (sol.t[end], false) | ||
| end |
There was a problem hiding this comment.
A fancier approach would be to use the integrator interface and keep stepping until a callback calls terminate! (assuming it will always be called eventually, but could with some non-zero probability be after tspan[2] -- I haven't checked the problem carefully enough to understand if that is the case here).
There was a problem hiding this comment.
I agree that this would be a cooler way. Maybe add a note about it for now? Now in principle that this is doable, but haven't actually seen it in practise.
| Finally, we can plot the distribution of activation times. For this, we will use the [StatsPlots.jl](https://docs.juliaplots.org/latest/generated/statsplots/) package's `density` function (essentially a smoothed histogram). The input to `density` is the activation times (which our output function has saved to `esol.u`). | ||
| ```@example activation_time_distribution_measurement | ||
| using StatsPlots | ||
| density(esol.u; label = "Activation time distribution", xlabel = "Activation time") | ||
| ``` |
There was a problem hiding this comment.
Why not just use a normal histogram with the option to have it normalize to a pdf? That shows the actual sample noise in the bins instead of smoothing everything out, and avoids the StatsPlots dependency.
There was a problem hiding this comment.
I generally have preferred to use density as I feel those plots are better and easier to work with. But if avoiding StatsPlots is priority I can plot a histogram and then link to the desnity function suggesting that readers should try it out as well?
…n_measurement.md Co-authored-by: Sam Isaacson <[email protected]>
Combines a decent couple of concepts from various parts of the docs, and also shows how to use
EnsembleProblemoutput functions.