Once I create ensemble systems, How could we save the solution in text or csv files ?
γ = 0.2647
D = 0.4811
μ= 34
u₀=34 #7.0
f(u,p,t) = -γ*(u - μ)
g(u,p,t) = sqrt(D)
dt = 1
tspan = (1,60)
prob = SDEProblem(f,g,u₀,tspan)
sol = solve(prob,EM(),dt=dt,save_noise=true)
ensembleprob = EnsembleProblem(prob)
sol = solve(ensembleprob,EnsembleThreads(),trajectories=100)
plotly()
plot(sol,linealpha=0.4,xlabel = “Time (years)”,ylabel = “Carbon Price (USD/tCO2)”)
Save it to a dataframe with:
data_save = DataFrame(sol)
CSV.write("Simulation.csv",data_save)
It seems that Ensemble problem can not simply saved via DataFrame()
This is the output:
ERROR: LoadError: type EnsembleSolution has no field prob
I am looking for a way to save each trajectory data so that I can do some post processings on it later.
Using BSON package also does not work, because in the solution field of each trajectory, I see that it changes solution to UInt8 format, while I expect an Array of Arrays.
Any suggestions?
If you check the output of ensemble trajectory, it is something like this:
size(sol) = (i,j,k)
where i is the number of variables in the variables in the ODE/SDE system, j is the number of steps in each trajectory, and k is the number of realizations/trajectories.
Having that in mind you can simply writhe output in this format:
for j = 1:trajectories
@inbounds solnew = vcat(j, sol.u[j].t',sol[:,:,j])
@inbounds open(out_file,"a+") do file
writedlm(file, solnew', ',')
end
end
I hope this work for you.
For hierarchical data like this, I highly suggest a format that directly supports hierarchy, but has a generic reader. This way, anyone can decipher the structure using the generic reader. HDF5 works great for this. You can download HDF5 view to see what’s inside (looks like files & folders).
My CMDimData.jl/MDDatasets.jl libraries provide multi-dimensional data structures and functions that can easily be written to/read from HDF5 files. An example of how you can transfer “Ensemble” results to these data structures is shown below.
Create ydataMD (of type: MDDatasets.DataRS) from “Ensemble” ydata array (and x values):
#Encapsulate multi-dimensional data using MDDatasets.jl structures:
n_trajectories = length(ydata)
ydataMD = fill(DataRS{DataF1}, PSweep("Trajectory", collect(1:n_trajectories))) do trajectory
sol_i = DataF1(x, ydata[trajectory]) #encapsulates (x,y) data in a single structure
return sol_i
end
If you are interested, but have trouble understanding how it works, I can provide further details.
Once ydataMD is created, you can easily post-process the data over all parameter sweeps without having to explicitly loop through each one. There is example code in the CMDimData.jl README file:
→ GitHub - ma-laforge/CMDimData.jl: Parametric analysis/visualization +continuous-f(x) interpolation
I recently made a few changes to CMDimData.jl, so the following code will only work once v0.4.0 is merged in Julia’s GeneralRegistry:
using MDDatasets #Data containers & post-processing operation
#Higher level plotting routines, supporting MDDataset objects:
using CMDimData, CMDimData.EasyPlot, CMDimData.EasyData
#Build some data like with your own "Ensemble" results:
x = collect(1:100)
ydata = []
for trajectory in 1:50
push!(ydata, rand(length(x)))
end
#Encapsulate multi-dimensional data using MDDatasets.jl structures:
n_trajectories = length(ydata)
ydataMD = fill(DataRS{DataF1}, PSweep("Trajectory", collect(1:n_trajectories))) do trajectory
sol_i = DataF1(x, ydata[trajectory]) #encapsulates (x,y) data in a single structure
return sol_i
end
#Read/write data to HDF5 file:
filepath = "results.hdf5"
datasetid = "EnsembleDataset" #Some string
@info("Writing $filepath...")
EasyData.openwriter(filepath) do w
write(w, ydataMD, datasetid)
end
@info("Reading back $filepath...")
r = EasyData.openreader(filepath)
ydataMD_fromfile = EasyData.readdata(r, datasetid)
close(r)
@show ydataMD_fromfile
As you can see, all data is written out with the write() and EasyData.readdata() commands (without having to manually loop through trajectories).
#Show plot (more control required than Plot.jl)
#------------------------------------------------------------------------------
plot = cons(:plot, title = "multi-dimensional data",
xyaxes = set(xscale=:lin, yscale=:lin),
labels = set(xaxis="X-Axis Label", yaxis="Y-Axis Label")
)
#Add all of ydataMD to the plot with a single statement:
push!(plot,
cons(:wfrm, ydataMD, label="ydataMD"),
)
#Display plot itself:
CMDimData.@includepkg EasyPlotInspect #Use InspectDR backend to plot (needs to run only once)
EasyPlot.displaygui(:InspectDR, plot) #Can use other backends
