_isinplace causes incorrect dispatch when both out of place and in place versions of a function is defined

[axla-io]

It seems that BifurcationKit._isinplacefails for some conditions, see MWE. A suggestion would be to change BifurcationKit._isinplace to SciMLBase.isinplace

MWE:

using BoundaryValueDiffEq
using NonlinearSolve
using BifurcationKit
using Test
tspan = (0.0, pi / 2)
function simplependulum!(du, u, p, t)
    θ = u[1]
    dθ = u[2]
    du[1] = dθ
    du[2] = -9.81 * sin(θ)
end
function bc!(residual, u, p, t)
    residual[1] = u[end ÷ 2][1] + pi / 2
    residual[2] = u[end][1] - pi / 2
end

bvp = BVProblem(simplependulum!, bc!, [pi / 2, pi / 2], tspan, dt = dt)
dt = 0.01

function get_nl_from_bvp(prob::BVDE.BVProblem, alg::BVDE.AbstractMIRK; dt=0.0,
    adaptive=true, kwargs...)
    cache = BVDE.init(prob, alg; dt = dt, adaptive = adaptive, kwargs...)
    @unpack y₀= cache
    return BVDE.construct_nlproblem(cache, BVDE.recursive_flatten(y₀))
end

nl_prob = get_nl_from_tpbvp(bvp, alg, dt = dt)

bf_prob = BifurcationProblem(nl_prob.f.f, nl_prob.u0, nl_prob.p, (@lens _[1]), J=nl_prob.f.jac) # BifurcationProblem is out of place but should be in place.

@test SciMLBase.isinplace(nl_prob.f.f, 3) == BifurcationKit._isinplace(nl_prob.f.f) # Expected to be true

[rveltz]

Hi

This part is undocumented on purpose. There is no inplace version of BifurcationProblem. The purpose of the BK.isinplace is to detect when a user passes an inplace version and build an out place version of it to simplify its life.

In the future id like to support inplace formulation but it is quite some work.

You can emulate this by doing

  const J = zeros(n,n)
MyJac!(x,p) = J .= 

And pass it to BifurcationProblem as J argument.

And

[axla-io]

Thanks, I see. Wouldn't it be nice though to change to SciMLBase.isinplace, so that the BifurcationProblem constructor could handle these kinds of cases too? Or is there a specific reason to use BifurcationKit._isinplace instead?

[rveltz]

The underscore _isinplace indicates it is an internal function to BK that the user should not use. You are right to think that, in the future, we will likely bind it to SciMLBase.isinplace when we support inplace problems. Just to give you a hint of the issue involved: KrylovKit (linear solver) wants out of place jacobian, IterativeSolvers accepts both,... Same for the EigenSolver.

I could perhaps use https://github.com/SciML/LinearSolve.jl to mitigate this (and in the mean time create a https://github.com/SciML/EigenSolve.jl package!) but BK also support types that are not AbtractArray (like KrylovKit). So you see, it involves decisions and quite some work.

Now, back to your problem, BifurcationProblem will be out of place no matter what you pass, but that is not really the performance bottleneck. It will be in the jacobian computation. Using a const J like I suggested above will mitigate this.

[axla-io]

Thank you for the clarification! I understand that there is a lot to go through in order to get the in place version working and that SciMLBase.isinplace will make sense later.

To solve my problem I did make an out of place Jacobian function.

[rveltz]

Let's keep it open.

To solve my problem I did make an out of place Jacobian function.

I am not sure I get your comment right. In case I did not: You can really use the inplace SciML jacobian using a gobal const for the return matrix and the computation of the jacobian will be allocation free. I hope this is what you do. Here is an example:

https://github.com/bifurcationkit/BifurcationKit.jl/blob/master/src/periodicorbit/PeriodicOrbitTrapeze.jl#L847

Here _J is re-written inplace and passed to BK.

[axla-io]

Oh I just did an allocating version of it... I never used a global const array before, because I was a bit scared of global variables, but I guess it's alright if it is a const. I tried your tip now and it gives a lot better performance than the allocating version. Thanks!

[rveltz]

By the way, the toy problem of this issue would be much better handled with orthogonal collocation (that is in BK)

[rveltz]

should we close this?

[axla-io]

Yeah, I think that we can!