ODE.jl
ODE.jl copied to clipboard
SciML
Change output from Array{Array{Float64,1},1} to Array{Float64,2}
When I run
t,x = ode45(f,x0,tspan)
I get:
typeof(x) # Array{Array{Float64,1},1}
Which is annoying because I want to, for example, plot the evolution for the first dynamic variable:
plot(x[:,1])
Which doesn't work... It just plots all the dynamic variables. Is there a workaround for this? Or an easy way to convert Array{Array{Float64,1},1} to Array{Float64,2}? Shouldn't ode45 be modified to return a 2D array?
There is a simple, but non-intuitive, way, to convert from a Vector{Vector{Float64}} to a Vector{Float64, 2}. It's as follows:
julia> v = Vector{Float64}[[1, 2], [3, 4], [5, 6]]
3-element Array{Array{Float64,1},1}:
[1.0,2.0]
[3.0,4.0]
[5.0,6.0]
julia> hcat(v...)'
3x2 Array{Float64,2}:
1.0 2.0
3.0 4.0
5.0 6.0
Here, hcat is a function that concatenates horizontally; then ' takes the transpose.
An alternative way to extract the first component from each is
xs = [x[1] for x in v]
ys = [x[2] for x in v]
plot(xs, ys)
Thanks a ton! This is very helpful.
How can we fix this? Would it be possible to add something like the following to Base:
convert(::Type{Array{T,2}}, v::Array{Array{T,1},1}) = hcat(v...)'
Or could we at least change the output of the ODE functions to give a Matrix{Float64} by default?
@ahwillia : Returning Vector{typeof(y0)} was actually a conscious choice for the API, discussed in #20. As Tomas explained recently in #78, the reason is that we also want to support "scalar-like types" which support some set of mathematical operations (see #37 for an example).
Last time I checked, using a list comprehension was faster. This is actually an unexported function in ODE.jl:
vcat_nosplat(y) = eltype(y[1])[el[1] for el in y] # Does vcat(y...) without the splatting
Anyway, I also find this a bit annoying. Could the 2D array be the standard output for 1D array IC input? And only the vector of vectors for scalar-like things?
Last, there is the scalar-like test https://github.com/JuliaLang/ODE.jl/blob/master/test/interface-tests.jl to look at.
Anyway, I also find this a bit annoying. Could the 2D array be the standard for 1D array IC input? And only the vector of vectors for scalar-like things?
To be honest I am not a big fan of such a change (but I am ready to be convinced otherwise). The current API might be (mildly) annoying, but it is intuitive and consistent: the output always is a vector of solutions of the type of the IC. You can use any of those solutions immediately as an IC. The price you have to pay is calling hcat(y...)' (or vcat_nosplat(y)) ...
How would you put the 1D arrays in the 2D output, row wise or column wise?
Yes, you're right. It's just that I almost always do a hcat/vcat, which somewhat is annoying too. I'd do a column per time step.
Without specific numbers, there might be some performance concern for using Vector{Vector{T}}. As a start, it allocates a lot more memories that are not necessarily contentious in memory. (assuming the temporary buffers are pre-allocated). It is also a bad idea to splat a huge array to vcat/hcat so if the API is to be kept this way, there should be an easier and better way to do the conversion (e.g. by exporting the vcat_nosplat)
How would you put the 1D arrays in the 2D output, row wise or column wise?
Column default makes sense and transposing should be much easier than flattening.
I don't really think consistency is an issue here. There's good reason that we have a nice type system and don't treat everything as a black box. A Vector or Array are distinct types from others and not returning a vector of them isn't a unreasonable special case.
Without specific numbers,
And now actual numbers. The benchmarks source can be found here and the results can be found in the same directory and is also pasted here,
ele_size = 5
create_vectors(ele_size,seq_len)
Time per evaluation: 264.36 μs [137.78 μs, 390.93 μs]
create_matrix(ele_size,seq_len)
Time per evaluation: 46.21 μs [34.85 μs, 57.56 μs]
transpose(res_matrix)
Time per evaluation: 65.40 μs [48.73 μs, 82.07 μs]
splat_vector1(res_vector)
Time per evaluation: 95.23 μs [80.84 μs, 109.63 μs]
splat_vector2(res_vector)
Time per evaluation: 126.31 μs [110.63 μs, 141.99 μs]
sum_vectors(res_vector)
Time per evaluation: 85.35 μs [84.95 μs, 85.75 μs]
sum_matrix(res_matrix)
Time per evaluation: 5.49 μs [5.48 μs, 5.50 μs]
ele_size = 100
create_vectors(ele_size,seq_len)
Time per evaluation: 901.62 μs [682.24 μs, 1.12 ms]
create_matrix(ele_size,seq_len)
Time per evaluation: 583.15 μs [506.91 μs, 659.40 μs]
transpose(res_matrix)
Time per evaluation: 1.18 ms [1.08 ms, 1.28 ms]
splat_vector1(res_vector)
Time per evaluation: 1.28 ms [1.17 ms, 1.39 ms]
splat_vector2(res_vector)
Time per evaluation: 4.34 ms [4.16 ms, 4.51 ms]
sum_vectors(res_vector)
Time per evaluation: 242.98 μs [231.56 μs, 254.41 μs]
sum_matrix(res_matrix)
Time per evaluation: 145.58 μs [134.14 μs, 157.03 μs]
The benchmark measures the effeciency for some basic operations including creating the buffers, converting them into different format and assessing elements (sum).
As clearly shown above, creating arrays of arrays is mush slower, Concatenating the vectors is commonly even more expensive than transposing (splat_vector1 is relatively effective because the memory is relatively continuous due to how it is created in this benchmark but it is still very expensive). There's no benchmark on converting from Matrix to Vector of Vector because you can just use SubArray or ArrayView to slice the matrix. Access the element in the buffer is also much more expensive when you have Vector of Vector. Not to mention that there's good reason we have a Matrix type and many libraries expect data in such format when getting a 2D array.
I am not saying that matrices are bad, it's just that Vector{typeofy0)} is a reasonable choice and gives a consistent API. I would also argue that the performance of the array operations you benchmarked is not critical if you are solving large systems of ODEs. However, I see that there are applications like solving of PDEs for which the Matrix output would be more convenient. If people think that this case sufficiently important, we could try to special-case Vector ICs and see how it goes ...
My point is, "consistency" is pretty much the only possible reason for the current behavior but it's not really a good reason because the other way around is consistent with the distinction between array and scalar every else.
For performance, you need a fairly big system for the overhead to not matter and it's very significant for small systems.
For using the result, the current one is almost certainly not the form you want. For actual ODE, at least for what I've dealt with, I usually want to see the evolution of each degrees of freedom. For PDE, as you said, it's handy to present the results as 2d arrays. Even if you are somehow more interested in each vector, simply constructing subarray or arrayview is easy and cheap.
The main argument for the current behavior is generality. I might want to integrate ODEs over some arbitrary Banach space, not just y::Vector. For example, I might want to apply Runge-Kutta to chebfun-like representations of functions. In general, there is no clean way to represent the result as a 2d array.
If you are doing arithmetics on a custom type, it's pefectly fine to return a Vector of that as the result. And that doesn't conflict with returning a Matrix (or a higher dimensional Array) when y0 is AbstractVector (AbstractArray) or at least Vector (or Array).
@yuyichao : Don't get me wrong, I don't insist on the current behavior. However, we went through some discussions (#20) to define the current API, so we shouldn't make fast decisions now.
My point is, "consistency" is pretty much the only possible reason for the current behavior but it's not really a good reason because the other way around is consistent with the distinction between array and scalar every else.
There is no logical reason why an ODE solver has to return a matrix either. Maybe this was introduced by MATLAB? Often solvers (like in Sundials) only return the final solution and one can access intermediate results only via some event system. The convenience really depends on the use case I suppose.
If you are saying that there is no Julia function/package that gives Vector{Vector{T}}, then we should definitely change the current behavior.
For using the result, the current one is almost certainly not the form you want.
That also depends on what you want to do. If you want to evaluate some function f(y) for each solution it is very easy to do with the current behavior (and this is probably not uncommon for PDEs).
so we shouldn't make fast decisions now.
Sure.
There is no logical reason why an ODE solver has to return a matrix either. Maybe this was introduced by MATLAB?
That would be a plus for least surprise. (Also for odeint in scipy I believe).
Often solvers (like in Sundials) only return the final solution and one can access intermediate results only via some event system. The convenience really depends on the use case I suppose.
I've never use Sundials but that interface sounds better. I also do this when I write these kind of functions to have a flexible interface to gather intermediate results, which could be a dummy one if you don't care about them or don't want to pay the overhead (or if they are simply too big to fit in the memory).
If you are saying that there is no Julia function/package that gives Vector{Vector{T}}, then we should definitely change the current behavior.
Well, Vector{Vector{T}} as a function returns the type Vector{Vector{T}} when called with a integer ;-p. More seriously though, I think most functions that expect 2-D data will want a Matrix or at least a AbstractMatrix. Accessing the history is in general harder when you have two indirections and cannot slice out a row easily.
That also depends on what you want to do. If you want to evaluate some function f(y) for each solution it is very easy to do with the current behavior (and this is probably not uncommon for PDEs).
Agree. The point is that it is equally easy to do with matrix. Just do res[:, i] if you can afford the copy or use sub(res, (:, i)) if you don't want to copy the data and before we turn array indexing into subarray yet. I think in general this is not true the other way around though, i.e. there are many things that are not as easy to do if you have a vector of vector instead of a matrix.
As another side note, I think the ode interface for scipy also allow flexible result gathering (via method overloading I believe).
This is above my paygrade, but fwiw I think the majority of users would appreciate having a Matrix as the default output, especially those (like myself) who are used to using odeint in scipy and MATLAB. I see the argument for handling more general use cases. Maybe multiple dispatch could allow for something more flexible, like @yuyichao suggests.
The plan was to use an iterator interface for the low level functions. Then
one could do for t,y in OdeXX(F,y0,tspan, args) ... end and have and
option for "event selection". This should offer the most flexibility - one
can store the result or just extract some component or whatever. The
current functions odexx would be written around the iterators.
On Do., 5. Nov. 2015 at 20:30, Alex Williams [email protected] wrote:
This is above my paygrade, but fwiw I think the majority of users would appreciate having a Matrix as the default output, especially those (like myself) who are used to using odeint in scipy and MATLAB. I see the argument for handling more general use cases. Maybe multiple dispatch could allow for something more flexible, like @yuyichao https://github.com/yuyichao suggests.
— Reply to this email directly or view it on GitHub https://github.com/JuliaLang/ODE.jl/issues/80#issuecomment-154164585.
Perhaps the problem is, really, that we're using Vector{typeof(y0)} as a shortcut for having a true "solution type" such as ODESolution{typeof(y0), <maybe more type params>}. It might be worth considering adding a type hierarchy for the output of our ODE solvers, to make it easier to optimize various cases.
For example, indexing into an ODESolution{Vector{T}} could be implemented so that when using linear indexing, it feels like a Vector{Vector{T}}, but multilinear indexing can make it feel like a Matrix{T}. The internal data structure can be whatever is the fastest, regardless.
This would probably result in quite a bit of extra code to maintain. I'm not sure it would be worth it.
Actually, some changes that made it into Julia 0.4 make it easier than ever to overload indexing; we only have to care about how indexing is handled for single elements, and then there are wrapper function in Base that handle indexing with ranges, colons etc.
I think something like this would suffice to start with:
abstract AbstractODESolution{T} <: AbstractVector{T}
immutable ODESolution{T} <: AbstractODESolution{T}
solution::Vector{T}
end
ODESolution{T}(::Type{T}) = ODESolution(Vector(T,0))
getindex(s::ODESolution, i::Integer) = getindex(s.solution, i)
push!{T}(s::ODESolution{T}, v::T) = push!(s.solution, v)
immutable VectorODESolution{T} <: AbstractODESolution{Vector{T}}
solution::Matrix{T} #
end
VectorODESolution{T}(::Type{T}) = VectorODESolution(Matrix(T,0,0))
ODESolution{T<:Vector}(::Type{T}) = VectorODESolution(Matrix(eltype(T),0,0)) # optional, but nice
getindex(s::VectorODESolution, i::Integer) = sub(s.solution, (:, i))
getindex(s::VectorODESolution, i::Integer, j::Integer) = getindex(s.solution, i, j)
push!{T}(s::VectorODESolution{T}, v::Vector{T} = hcat(s.solution, v) # is there a way to avoid this copy?
As you notice, the hcat at the end actually copies the solution each time a vector is added, which is obviously going to be bad for performance. This isn't unique for this approach, though; if we decide to just return a Matrix directly, how would you build it iteratively while solving the ODE?
This isn't unique for this approach, though; if we decide to just return a Matrix directly, how would you build it iteratively while solving the ODE?
You can append to a vector and reshape if the size isn't known beforehand.
Ah! That makes sense. That's also good because it can leave the internals of the solvers untouched, using Vector{typeof(y0)}s until we reach the end of the run, and then wrap it in a subtype of AbstractODESolution. Both push! methods above can be discarded.
@tlycken +1 on the idea of AbstractODESolution. Having solutions types is useful not just for the issue here but could also reflect the fact that ODEs can have "special" solutions in certain cases (for example for particular initial or final conditions it might be easy to find the solution).
But with respect to this issue, there could be a benefit in dropping the AbstractVector constraint to allow for more general behavior. For a first example, in cases where an exact solution is known as a function of time it might be more useful to have a callable type that can return the known solution at arbitrary times. Building on this first example: for a scalar ode y' = 0 + g(t) where g is a polynomial the solution will be polynomial an could in fact be represented by a polynomial type if it were advantageous to do so.
As far as I can tell removing the constraint doesn't cost anything, since it can always be put back in later (e.g. with some kind of AbstractODESolutionVector<:AbstractODESolution).
I'd also like to +1 @stevengj comment about allowing for generality in the type of the state variable. If his idea about chebfun-like representations sounds too exotic to some, I'd like to point out that it would also be desirable to accommodate the case where the state variable is a matrix (i.e. matrix differential equation) or higher order array, the applications of which are quite common. Having to pack/unpack matrices into vectors to set solutions to matrix odes was one of my number one pet peeves with fortran and MATLAB. Julia could lead to way to a better future for this area:)