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

Question

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

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

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

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 будут напечатаны без ошибок преобразования.

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

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