Histogram from a set of values from EnsembleProblem

New to Julia
1

Dear All,
I have a system of six SDEs I solved with

prob_sde = SDEProblem(qdelta,qsigma,u0,(0.0,600.0),noise=W)
ensembleprob = EnsembleProblem(prob_sde,prob_func=prob_func)
sol = solve(ensembleprob,SRIW1(),EnsembleDistributed(),trajectories=50)
plotly()
plot(sol,vars=(0,6),linealpha=0.4)

which seems to work well, unfortunately I need the histogram of u[6], but I merely don’t know how to extract from ‘sol’ all the values from the 50 trajectories so that I can simply call the histogram ploting tool… I tried convert(Array,sol) which failed, then v=hcat(sol.u[6,:]) which failed too…

sorry again for this naive question
thank you !
Frederic

2

Each trajectory has a different number of points due to adaptivity. You either need to adjust your thinking about the timeseries to that or just make it save at even spots, like saveat=1 so that each trajectory has the same saving spots.

3

thank you for the help ! the saveat option will give me values at a given time for all trajectories, right ? I guess it is not adequate for my case.
let suppose I need the histogram of only the nth trajectory, how do I extract the values then as hcat refuse to do so ?
sorry for my basics questions !

4

v=reduce(hcat,[sol[i].u[6,:] for i in 1:length(sol)]). Remember sol is an ensemble of trajectories, so you want to grab the timeseries of the 6th variable from each trajectory, and then hcat them all together (via reduce(hcat,x)).

5

Thanks Chris ! the outuput of v is 1×50 Array{Array{Float64,1},2} so when I just do histogram(v[1,:]) result is … strange.

So I am suppose to have 50 timesseries of a lot of points (600seconds), but obviously the reduce reduced them all :slight_smile: .

what I don’t get is when I ask to plot(sol[1],vars=(0,6)), it works well

Capture d’écran 2020-04-07 à 12.18.19
Capture d’écran 2020-04-07 à 12.18.191194×778 198 KB

I tried to read and understand the manual page on MultiDimensional arrays but is doesn’t help me much…

6

hello
I’m coming back to this problem. As far as I understood, the v=reduce(hcat,[sol[i].u[6,:] for i in 1:length(sol)]) output very low amount of data, so I think it is more u[:,6] which correspond to plot(sol[1],vars=(0,6)), but when I do so Julia does not agree :

BoundsError: attempt to access 13796-element Array{Array{Float64,1},1} at index [Base.Slice(Base.OneTo(13796)), 6]

So how is it possible to extract, for the solution i, all the values of the 6th variable and use it in the GR histogram function ?

Sorry for such basic question but at that time I’m stuck !

7

sol[1,6,:]

8

thanks a lot ! in fact sol[1,:,6] in my case no ? else I have only 20 elements. then histogram works well.

9

Oh, sol[6,:,1] my bad.

10

then I am a bit lost again, what are the i,j,k in sol[i,j,k] ?
initially I was hoping for
i → trajectory
j → time
k → variable
to comply with sol[i],vars=(j,k)

11

No, it’s column order, it’s

i → variable
j → time
k → trajectory

s[i,j] is variable and time, and then s[i,j,k] is variable time trajectory. That’s the only thing that makes sense for column-ordering since trajectories are stored in memory together, so the 3rd axis is trajectory. Arrays in the language are made to be column order so if you orient it that way it’s clear.

12

ok understood, anyway if I call  ht = sol[6,:,15] it output BoundsError: attempt to access 14223-element Array{Array{Float64,1},1} at index [14224] I am supposed to have 20 trajectories,
sol = solve(ensembleprob,SRIW1(),EnsembleDistributed(),trajectories=20);
but only k=1 and 2 works… where am I wrong ?

13

Try sol = solve(ensembleprob,SRIW1(),EnsembleDistributed(),trajectories=20,saveat=0.1); or something like that: the issue is raggedly-sized trajectories since adaptive timestepping allows each different trajectory to have different amounts of outputs (which is why hcating failed)

14

In that case it works as every trajectory has the same size. Thanks a lot !

What I don’t get is why Julia doesn’t want to output one array even if its size is different from other arrays… why is it compared to others ? As it exists it can be extracted no ? worst process ever might be element by element and rebuilding an array that way …

15

It’s not Julia, it’s my definition of indexing on a VectorOfArray. All of the indexing stuff somewhat assumes size(x) makes sense, which it doesn’t when the outputs are ragged. It would be good to get an issue on this.

16

Track it here: Indexing ragged size EnsembleSolutions · Issue #597 · SciML/DifferentialEquations.jl · GitHub