Юлия Дифференциальные уравнения: не удается решить проблему MonteCarloProblem, когда работает «ручное» решение (ошибка преобразования)

Question

Я работаю над немного большим проектом, в котором я интенсивно использую типы MonteCarloProblem DifrentialEquations.jl.Обычно все работает нормально, но я столкнулся со странным случаем, в котором я не могу разобраться, потому что ручная итерация по каждой отдельной проблеме работает просто отлично.

Я столкнулся с этим, когда опробовал некоторые системы, в которых решения могут быстро вырасти до бесконечности, и я наивно использовал интегратор с постоянной шириной шага, такой как Эйлер, чтобы решения могли содержать много инф.

Небольшой пример - хорошо известная 2-мерная нормальная форма бифуркации седлового узла с произвольной выборкой параметра и начальных условий.Таким образом, мы получим много инф в этом случае.Я знаю, что это не кажется особенно полезным, но это простой минимальный пример, который я мог бы придумать.

using DifferentialEquations
using Distributions
using StatsBase

struct saddle_node_parameters
   mu 
end

function saddle_node(du,u,p::saddle_node_parameters,t)
   du[1] = p.mu - u[1]*u[1]
   du[2] = - u[2]
end

N = 50
ICs = rand(Uniform(-2,2),2,N)
pars = rand(Uniform(-1,1),N)
odep = ODEProblem(saddle_node, ICs[:,1], (0.,100.), saddle_node_parameters(0.5))
new_prob_func = (prob,i,repeat) -> ODEProblem(saddle_node, ICs[:,i],(0.,100.), saddle_node_parameters(pars[i]))
eval_func = (sol,i) -> (mean_and_std(sol[1,:]),false)
mcp = MonteCarloProblem(odep, prob_func=new_prob_func, output_func=eval_func)

Решение этой проблемы с помощью алгоритма по умолчанию (или автоматического определения алгоритма) работает (конечно, вы получите много предупреждений по мере того, как решения уходят в бесконечность, но без прерывающих ошибок).

Решение с постоянной шириной шага Эйлер не работает:

solve(mcp, alg=Euler(), dt=0.1, num_monte=N)

приводит к

MethodError: Cannot `convert` an object of type Nothing to an object of type Array{Array{Float64,1},1}

...

Stacktrace:
 [1] (::getfield(Base, Symbol("##683#685")))(::Task) at ./asyncmap.jl:178
 [2] foreach(::getfield(Base, Symbol("##683#685")), ::Array{Any,1}) at ./abstractarray.jl:1835
 [3] maptwice(::Function, ::Channel{Any}, ::Array{Any,1}, ::UnitRange{Int64}) at ./asyncmap.jl:178
 [4] #async_usemap#668 at ./asyncmap.jl:154 [inlined]
 [5] #async_usemap at ./none:0 [inlined]
 [6] #asyncmap#667 at ./asyncmap.jl:81 [inlined]
 [7] #asyncmap at ./none:0 [inlined]
 [8] #pmap#215(::Bool, ::Int64, ::Nothing, ::Array{Any,1}, ::Nothing, ::Function, ::Function, ::Distributed.CachingPool, ::UnitRange{Int64}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.0/Distributed/src/pmap.jl:126
 [9] #pmap at ./none:0 [inlined]
 [10] solve_batch(::MonteCarloProblem{ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,saddle_node_parameters,ODEFunction{true,typeof(saddle_node),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},getfield(Main, Symbol("##27#28")),getfield(Main, Symbol("##29#30")),getfield(DiffEqBase, Symbol("##378#384")),Array{Any,1}}, ::Nothing, ::Symbol, ::UnitRange{Int64}, ::Int64, ::Pair{Symbol,Any}, ::Vararg{Pair{Symbol,Any},N} where N) at /home/max/.julia/packages/DiffEqMonteCarlo/c9ztK/src/solve.jl:47
 [11] macro expansion at ./logging.jl:316 [inlined]
 [12] #solve#1(::Int64, ::Int64, ::Int64, ::Symbol, ::Base.Iterators.Pairs{Symbol,Any,Tuple{Symbol,Symbol},NamedTuple{(:alg, :dt),Tuple{Euler,Float64}}}, ::Function, ::MonteCarloProblem{ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,saddle_node_parameters,ODEFunction{true,typeof(saddle_node),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},getfield(Main, Symbol("##27#28")),getfield(Main, Symbol("##29#30")),getfield(DiffEqBase, Symbol("##378#384")),Array{Any,1}}, ::Nothing, ::Type{Val{true}}) at /home/max/.julia/packages/DiffEqMonteCarlo/c9ztK/src/solve.jl:21
 [13] (::getfield(DiffEqBase, Symbol("#kw##solve")))(::NamedTuple{(:alg, :dt, :num_monte),Tuple{Euler,Float64,Int64}}, ::typeof(solve), ::MonteCarloProblem{ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,saddle_node_parameters,ODEFunction{true,typeof(saddle_node),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},getfield(Main, Symbol("##27#28")),getfield(Main, Symbol("##29#30")),getfield(DiffEqBase, Symbol("##378#384")),Array{Any,1}}, ::Nothing, ::Type{Val{true}}) at ./none:0 (repeats 3 times)
 [14] top-level scope at In[25]:1

, даже если он работает, когда я просто вручную перебираю все проблемыкак я думаю, MonteCarloProblem делает это тоже:

for i=1:N 
   soli=solve(mcp.prob_func(odep, i, false),alg=Euler(), dt=0.1)
   println(eval_func(soli,i))
end

приводит к тому, что все наборы средств и sds печатаются просто отлично, без ошибок преобразования.

Я знаю, что этот пример может показаться не очень полезным, но для моего реального проекта я хотел бы иметь последовательное поведение, даже когда решения вырастают до бесконечности.

Заранее спасибо за любую помощь.