Skip to content
New issue

Have a question about this project? # for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “#”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? # to your account

histogram() function for result of ensemble simulation #657

Open
ignace-computing opened this issue Mar 8, 2023 · 5 comments
Open

histogram() function for result of ensemble simulation #657

ignace-computing opened this issue Mar 8, 2023 · 5 comments

Comments

@ignace-computing
Copy link
Contributor

Hello,
Sorry if this is not the right place to ask.
I am willing to contribute if the answer to the following question is negative. (In that case, please indicate some first steps.)
Let's say
sim = solve(prob,alg,ensemblealg,kwargs...)
Is there a function to plot a histogram of sim at a given simulation time, available?
Thanks and cheers,

@ChrisRackauckas
Copy link
Member

There is not a function for this right now.

@ignace-computing
Copy link
Contributor Author

Ok, I will make one.

@ignace-computing
Copy link
Contributor Author

Where are the Plot recipes for DifferentialEquations implemented?
I am having a hard time finding the location of this source code.

@ignace-computing
Copy link
Contributor Author

Ok, sorry, found it: SciMLBase/src/solutions/solution_interface.jl, for instance.

@ignace-computing
Copy link
Contributor Author

A first attempt is added below, I am still learning...

# here is a function that plots a histogram of all components of sim
# TODO: add meaningful labels
using DifferentialEquations.EnsembleAnalysis

function histogram_timepoint(sol::EnsembleSolution, t)
	return histogram(componentwise_vectors_timepoint(sol, t))
end

# here is a recipe to be used if you want to see a marginal distribution
# TODO: add labels 
using RecipesBase
@recipe function f(sol::EnsembleSolution, t; idx=nothing)
	componentwise_vectors_timepoint(sol, t)[idx]
end

## MWE 
# example copied from 
# https://docs.sciml.ai/DiffEqDocs/stable/features/ensemble/#Example-4:-Using-the-Analysis-Tools
function f(du,u,p,t)
  for i = 1:length(u)
    du[i] = 1.01*u[i]
  end
end
function σ(du,u,p,t)
  for i in 1:length(u)
    du[i] = .87*u[i]
  end
end
using DifferentialEquations
prob = SDEProblem(f,σ,ones(4,2)/2,(0.0,1.0)) #prob_sde_2Dlinear

prob2 = EnsembleProblem(prob)
sim = solve(prob2,SRIW1(),dt=1//2^(3),trajectories=100,adaptive=false)

using Plots 
p1 = histogram_timepoint(sim, 0.1)
 
p2 = histogram(sim, 0.1, idx=(1)) # using the above-defined recipe

If I add the above functions to SciMLBase, is this 'enough' for a pull request (after I will have added the labels using SciMLBase's add_labels! function, for instance)?
Any feedback is welcome, of course.

# for free to join this conversation on GitHub. Already have an account? # to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants