No method matching zero(::Type{Array{Float64,1}})

New to Julia
1

I am trying to implement a Household epidemiology model. Below, I provide a minimal working example to demonstrate the error (of course, the example is no longer meant to model the epidemic):

using DifferentialEquations
using SimpleDiffEq
# Use commas to avoid printing in Atom when writing code
nb_nodes = 4
S0 = Array{Float64,1}(undef, 4);
fill!(S0, 4);

function households!(du,u,p,t)
	(S) = u
	(βG) = p
	infections = @. βG * sum(I) + βG * I
	du[1] = -infections .* S
	nothing
end;  # no output due to semi-colon

u0 = [S0]
params = [0.5]
prob_ode = ODEProblem(dbug!,u0,tspan,params);
sol_ode = solve(prob_ode);

When executing solve(prob_ode), I get the following error, which I do not know how to debug. Can anybody help? Here is the error:

image
image722×277 42.4 KB

2

The error message is saying that it doesn’t know how to calculate a zero value of the type Vector{Float64}.

But I am unable to reproduce that error. You seem to be missing the definition of the dbug! function.

Also note that the two lines for defining S0 can be reduced to S0 = fill(4.0, 4).

Perhaps you want u0 = S0 instead of u0 = [S0]?

1 Like
3

Sorry about that. The undefined function was defined in another Atom window. I exited and restarted Atom. Try this code, which reproduces the error. I have a feeling that I cannot define systems of matrix equations with DifferentialEquations.jl .

using DifferentialEquations
using SimpleDiffEq
S0 = Array{Float64,1}(undef, 4);
fill!(S0, 4);

function households!(du,u,p,t)
	(S) = u
	(βG) = p
	infections = βG * sum(I) .+ βG .* S
	du[1] = -infections .* S
	nothing
end;  

tspan = (0.,5.)
u0 = [S0]
params = [0.5]
prob_ode = ODEProblem(households!,u0,tspan,params);
sol_ode = solve(prob_ode, Tsit5());```

The error is as before;

image
image722×277 42.4 KB

Here is the problem as I see it: my variables are described as a list of arrays, and therefore, so is the u0 vector. DifferentialEquations.jl probably uses zeros() to initialize u0, and zeros cannot be applied to an array. Bottom line, it would seem that I cannot define systems of arrays. If this is true, it would be nice to extend DifferentialEquations.jl to handle this case. In the meantime, what is the solution? Is it to create one large matrix system, i.e. combine all my equations into one large equation?

4

If they are an array of small fixed-size arrays, as in the 4-element arrays in your example above, consider using an array of SVector (from StaticArrays.jl) instead. That way, the length of the arrays are known to the compiler and zero will work. And, incidentally, it should be much faster.

That being said, it might be possible to fix DifferentialEquations.jl to work with your case. It’s hard to tell, because you posted a screenshot that doesn’t show the full backtrace, so it’s not clear precisely where the error comes from. Don’t post screenshots (they aren’t searchable, copyable, …).

3 Likes
5

Sorry about that. Here is my example again. The reason that my arrays are size 4 is because this is a MWE. So StaticArrays are not really appropriate in this case. Notice that u0 is formed from two arrays since this is more representative of the actual problem I am working with. In the meantime, I have concatenated the two arrays into one long array and am using views.

using DifferentialEquations
using SimpleDiffEq

# Use commas to avoid printing in Atom when writing code
S0 = Array{Float64,1}(undef, 4);
I0 = Array{Float64,1}(undef, 4);
fill!(S0, 4);
fill!(I0, 4);

function households!(du,u,p,t)
	(S,I) = u
	(βG) = p
	infections = βG * sum(S) .+ βG .* I
	du[1] = -infections .* S
	nothing
end;  # no output due to semi-colon

tspan = (0.,5.)
u0 = [S0, I0]
params = [0.5]
prob_ode = ODEProblem(households!,u0,tspan,params);
sol_ode = solve(prob_ode, Tsit5());

Here is the complete error trace:

MethodError: no method matching zero(::Type{Array{Float64,1}})
Closest candidates are:
  zero(!Matched::Type{Missing}) at missing.jl:103
  zero(!Matched::Type{LibGit2.GitHash}) at /Users/sabae/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/LibGit2/src/oid.jl:220
  zero(!Matched::Type{Pkg.Resolve.VersionWeight}) at /Users/sabae/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/Pkg/src/Resolve/versionweights.jl:15
  ...
zero(::Array{Array{Float64,1},1}) at abstractarray.jl:902
alg_cache(::Tsit5, ::Array{Array{Float64,1},1}, ::Array{Array{Float64,1},1}, ::Type{T} where T, ::Type{T} where T, ::Type{T} where T, ::Array{Array{Float64,1},1}, ::Array{Array{Float64,1},1}, ::ODEFunction{true,typeof(households!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing}, ::Float64, ::Float64, ::Float64, ::Array{Float64,1}, ::Bool, ::Val{true}) at low_order_rk_caches.jl:349
__init(::ODEProblem{Array{Array{Float64,1},1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(households!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::Tsit5, ::Array{Array{Array{Float64,1},1},1}, ::Array{Float64,1}, ::Array{Any,1}, ::Type{Val{true}}; saveat::Array{Float64,1}, tstops::Array{Float64,1}, d_discontinuities::Array{Float64,1}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Bool, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Nothing, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Nothing, reltol::Nothing, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, qoldinit::Rational{Int64}, fullnormalize::Bool, failfactor::Int64, beta1::Nothing, beta2::Nothing, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at solve.jl:270
__init(::ODEProblem{Array{Array{Float64,1},1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(households!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::Tsit5, ::Array{Array{Array{Float64,1},1},1}, ::Array{Float64,1}, ::Array{Any,1}, ::Type{Val{true}}) at solve.jl:67
__init(::ODEProblem{Array{Array{Float64,1},1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(households!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::Tsit5, ::Array{Array{Array{Float64,1},1},1}, ::Array{Float64,1}, ::Array{Any,1}) at solve.jl:67
__init(::ODEProblem{Array{Array{Float64,1},1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(households!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::Tsit5, ::Array{Array{Array{Float64,1},1},1}, ::Array{Float64,1}) at solve.jl:67
__init(::ODEProblem{Array{Array{Float64,1},1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(households!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::Tsit5, ::Array{Array{Array{Float64,1},1},1}) at solve.jl:67
__init(::ODEProblem{Array{Array{Float64,1},1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(households!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::Tsit5) at solve.jl:67
__solve(::ODEProblem{Array{Array{Float64,1},1},Tuple{Float64,Float64},true,Array{Float...
6

Arrays of arrays aren’t going to work so well because they don’t have a well-defined norm. You can use VectorOfArray from RecursiveArrayTools as the initial condition, and I think that should work because it has a norm and it has broadcast (don’t think I’ve tried it though).

u0 as a matrix might be easier here though.

2 Likes
7

u0 was the easy part though. I used

u0 = vcat(S0, I0, R0)
8

you mean hcat to make it a matrix, and then S = @view u[:,1]?

9

Here is my code, so you can see. I used vcat.


# Simulate Nh households of constant size n
# Assume a constant global transmission βG,
# and a local transmission βL

using DifferentialEquations
using SimpleDiffEq
using DataFrames
using StatsPlots
using BenchmarkTools
10

What do you mean? It is easy to define a norm for arrays of arrays, and the built-in LinearAlgebra.norm function does so:

julia> norm([[3,4,5], [6,7,8]])
14.106735979665883

(The fact that broadcast is not recursive may be a problem for you, though.)

11

Oh, I didn’t know that worked. Last time I tried this (admittedly years ago) there was issues with things recursing, and IIRC it was broadcast and norms. Indeed, defining a norm is simple, but it looks like that should work by default now which is nice.

But yes, that’s the second issue I remember, and this one is much more difficult. When I looked into this I opened an issue in Julia’s Base:

For other unrelated reasons though, I think @YingboMa is going to add recursive broadcast to RecursiveArrayTools.VectorOfArray soon (and it’ll already work if your internal vectors are all the same size), so u0 as a RecursiveArrayTools.VectorOfArray is a way around this until such a feature exists.

1 Like
12

I think you meant to post more.

13

I do not know what happened.

Here is the code:


# Simulate Nh households of constant size n
# Assume a constant global transmission βG,
# and a local transmission βL

using DifferentialEquations
using SimpleDiffEq
using DataFrames
using StatsPlots
using BenchmarkTools
14

EDIT: Didn’t notice the VectorOfArray suggestion above. I forgot about what I said below. You should probably just use that because you don’t have to change the function you’ve already written.

If you want to pass nested arrays through the ODE solver, you can use either `ArrayPartition`s--which are already exported in DifferentialEquations.jl--or `ComponentArray`s--which are from [ComponentArrays.jl](https://github.com/jonniedie/ComponentArrays.jl). `ArrayPartition`s are the closest to matching the syntax you've already written, but `ComponentArray`s have the added benefit of accessing subarrays through a named index. Here is how you'd do your above example in both:

ArrayPartitions:

using DifferentialEquations

# Use commas to avoid printing in Atom when writing code
S0 = Array{Float64,1}(undef, 4);
I0 = Array{Float64,1}(undef, 4);
fill!(S0, 4);
fill!(I0, 4);

function households!(du,u,p,t)
	S, I = u.x
	βG = p
	infections = βG * sum(S) .+ βG .* I
	du.x[1] .= -infections .* S
	nothing
end;  # no output due to semi-colon

tspan = (0.,5.)
u0 = ArrayPartition(S0, I0)
params = [0.5]
prob_ode = ODEProblem(households!,u0,tspan,params);
sol_ode = solve(prob_ode, Tsit5());

ComponentArrays:

using DifferentialEquations
using ComponentArrays

# Use commas to avoid printing in Atom when writing code
S0 = Array{Float64,1}(undef, 4);
I0 = Array{Float64,1}(undef, 4);
fill!(S0, 4);
fill!(I0, 4);

function households!(du,u,p,t)
	S, I = u.S, u.I
	βG = p
	infections = βG * sum(S) .+ βG .* I
	du.S = -infections .* S
	nothing
end;  # no output due to semi-colon

tspan = (0.,5.)
u0 = ComponentArray(S=S0, I=I0)
params = [0.5]
prob_ode = ODEProblem(households!,u0,tspan,params);
sol_ode = solve(prob_ode, Tsit5());

Full disclosure though: I wrote the ComponentArrays.jl package. I’m obviously a little biased here.

15

ComponentArrays.jl is a good solution as well. I think a lot of it is an X-Y problem: the really isn’t how to make arrays of arrays work, it’s how to get a nice abstraction for things which naturally have clumps with names, and ComponentArrays.jl, ArrayPartitions, LabelledArrays.jl, MultiScaleArrays.jl, etc. are all nice ways to do this.

2 Likes