Fix resize_J_W! for Vector{Matrix{Float64}}

[utkarsh530]

Attempts to Fix #1436.

MWE:

f = function (du,u,p,t)
  for i in 1:length(u)
    du[i] = (0.3/length(u))*u[i]
  end
end

condition = function (u,t,integrator)
  1-maximum(u)
end

affect! = function (integrator)
  u = integrator.u
  resize!(integrator,length(u)+1)
  maxidx = findmax(u)[2]
  Θ = rand()/5 + 0.25
  u[maxidx] = Θ
  u[end] = 1-Θ
  nothing
end

callback = ContinuousCallback(condition,affect!)

u0 = [0.2]
tspan = (0.0,10.0)
prob = ODEProblem(f,u0,tspan)

sol = solve(prob, ImplicitHairerWannerExtrapolation(),callback=callback)

Previously:

ERROR: MethodError: no method matching Vector{Matrix{Float64}}(::Matrix{Matrix{Float64}})
The closest candidates are:
  Array{T, N}(::AbstractArray{S, N}) where {T, N, S} at array.jl:540
  Vector{T}() where T at boot.jl:467
  Vector{T}(::Core.Compiler.AbstractRange{T}) where T at range.jl:1042

After:

┌ Warning: Interrupted. Larger maxiters is needed.
└ @ SciMLBase ~/.julia/packages/SciMLBase/cU5k7/src/integrator_interface.jl:331
retcode: MaxIters
Interpolation: 3rd order Hermite
t: 196073-element Vector{Float64}:
 0.0
 0.4237279064034065
 4.661006970437471
 5.385928365605549
 5.385928365605549
 5.38592836560555
 ⋮
 5.385928365779689
 5.38592836577969
 5.385928365779691
 5.385928365779692
 5.385928365779693
u: 196073-element Vector{Vector{Float64}}:
 [0.2]
 [0.2271102854460072]
 [0.8096636848795802]
 [0.9999999999999999]
 [0.630100191321227]
 [0.6301001913212272]
 ⋮
 [0.6301001913212272]
 [0.6301001913212272]
 [0.6301001913212272]
 [0.6301001913212272]
 [0.6301001913212272]

The problem we currently have with Extrapolation Methods and QNDF is the same; resizing of integrator.u is not happening.

See in QNDF:

 Warning: Interrupted. Larger maxiters is needed.
└ @ SciMLBase ~/.julia/packages/SciMLBase/cU5k7/src/integrator_interface.jl:331
retcode: MaxIters
Interpolation: 3rd order Hermite
t: 27-element Vector{Float64}:
 0.0
 3.350000000000003e-5
 6.700000000000006e-5
 0.00010050000000000008
 0.00043550000000000034
 0.0007705000000000006
 ⋮
 3.7028119459535453
 4.44080242793806
 5.178792909922574
 5.357995644079828
 5.357995644079828
u: 27-element Vector{Vector{Float64}}:
 [0.2]
 [0.20000169621691732]
 [0.20000365726070232]
 [0.20000565966574863]
 [0.20002539992745905]
 [0.2000455165235902]
 ⋮
 [0.6083341846429979]
 [0.7592640782476203]
 [0.9476395095023401]
 [0.9999999999999998]
 [0.5845412151714129]

So yeah, we need to address it.

[ChrisRackauckas]

this would alias all of them to the same W matrix. It needs to be an array of W's.

[ChrisRackauckas]

Needs tests.

[ChrisRackauckas]

Ping me to merge when this passes.

[utkarsh530]

I will be adding tests shortly, but with the given MWE (which is in ode_cache_tests.jl), why resizing of integrator.u does not happen in the cases of QNDF and Extrapolation Methods?

[ChrisRackauckas]

I don't get the question.

[utkarsh530]

Okay yeah, so in tests of ode_cache_tests.jl (prob is the same as defined above):

julia> sol = solve(prob,KenCarp4(),callback=callback)
julia> sol = solve(prob,KenCarp4(),callback=callback)
retcode: Success
Interpolation: 3rd order Hermite
t: 11-element Vector{Float64}:
  0.0
  0.07607837045701699
  0.8368620750271869
  2.179626217247377
  3.6441612528294822
  5.365296320151268
  5.365296320151268
  7.361064902152912
  7.588590959526431
  7.588590959526431
 10.0
u: 11-element Vector{Vector{Float64}}:
 [0.2]
 [0.20461719206588685]
 [0.2570769529714509]
 [0.3845969160080552]
 [0.5967754747467836]
 [0.9999999999999998]
 [0.28357714068834405, 0.716422859311656]
 [0.3825454441596969, 0.9664541375101127]
 [0.39582369127753253, 0.9999999999999998]
 [0.39582369127753253, 0.36257658405033333, 0.6374234159496667]
 [0.5037646054406076, 0.46145102941312277, 0.8112484490205522]

julia> sol = solve(prob,ImplicitHairerWannerExtrapolation(),callback=callback)
┌ Warning: Interrupted. Larger maxiters is needed.
└ @ SciMLBase ~/.julia/packages/SciMLBase/cU5k7/src/integrator_interface.jl:331
retcode: MaxIters
Interpolation: 3rd order Hermite
t: 196073-element Vector{Float64}:
 0.0
 0.4237279064034065
 4.661006970437471
 5.385928365605549
 5.385928365605549
 5.38592836560555
 ⋮
 5.385928365779689
 5.38592836577969
 5.385928365779691
 5.385928365779692
 5.385928365779693
u: 196073-element Vector{Vector{Float64}}:
 [0.2]
 [0.2271102854460072]
 [0.8096636848795802]
 [0.9999999999999999]
 [0.6428414471456899]
 [0.6428414471456898]
 ⋮
 [0.6428414471456899]
 [0.6428414471456899]
 [0.6428414471456899]
 [0.6428414471456899]
 [0.6428414471456899]

julia> sol = solve(prob,QNDF(),callback=callback)
┌ Warning: Interrupted. Larger maxiters is needed.
└ @ SciMLBase ~/.julia/packages/SciMLBase/cU5k7/src/integrator_interface.jl:331
retcode: MaxIters
Interpolation: 3rd order Hermite
t: 27-element Vector{Float64}:
 0.0
 3.350000000000003e-5
 6.700000000000006e-5
 0.00010050000000000008
 0.00043550000000000034
 0.0007705000000000006
 ⋮
 3.7028119459535453
 4.44080242793806
 5.178792909922574
 5.357995644079828
 5.357995644079828
u: 27-element Vector{Vector{Float64}}:
 [0.2]
 [0.20000169621691732]
 [0.20000365726070232]
 [0.20000565966574863]
 [0.20002539992745905]
 [0.2000455165235902]
 ⋮
 [0.6083341846429979]
 [0.7592640782476203]
 [0.9476395095023401]
 [0.9999999999999998]
 [0.556273938245464]

Both of these algorithms fail to solve contrary in algo like KenCarp4. Why is that happening?

[ChrisRackauckas]

Probably because W and J are aliasing instead of creating a new array of unaliased matrices.

[utkarsh530]

Where should be the tests added for this?

[ChrisRackauckas]

To the cache tests. You should be able to flip the currently failing tests.