Change output from Array{Array{Float64,1},1} to Array{Float64,2}

[ahwillia]

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?

[dpsanders]

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)

[ahwillia]

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?

[acroy]

@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).

[mauro3]

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.

[acroy]

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?

[mauro3]

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.

[yuyichao]

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.

[yuyichao]

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.

[acroy]

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 ...

[yuyichao]

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.

[stevengj]

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.

[yuyichao]

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).

[acroy]

@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).

[yuyichao]

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.

[yuyichao]

As another side note, I think the ode interface for scipy also allow flexible result gathering (via method overloading I believe).