internalnorm does not work with inplace methods

[Datseris]

MWE:

using OrdinaryDiffEq

function eom(du, u, p, t)
    @inbounds for i in 1:4
        du[1, i] = u[3, i]
        du[2, i] = u[4, i]
        du[3, i] = -u[1, i] - 2u[1, i]*u[2, i]
        du[4, i] = -u[2, i] - (u[1, i]^2 - u[2, i]^2)
    end
    return nothing
end

prob = ODEProblem{true}(eom, rand(4,4), (0.0, Inf))

# works
integ = init(prob, alg = Tsit5(), abstol = 1e-6, reltol = 1e-6)
# does not work:
integ = init(prob, alg = Tsit5(),
abstol = 1e-6, reltol = 1e-6, internalnorm = u -> norm(u[:, 1]))

it also doesn't work If I make the equations of motion a vector (it has nothing to do with them being in matrix form)

[devmotion]

IMO the problem is that internalnorm is used in places where we should rather just use abs, such as https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/initdt.jl#L11 and https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/initdt.jl#L125 (which should be broadcasted BTW to be consistent with the in-place computations) and also in https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/misc_utils.jl#L59 (and similar lines in the same file). For the default ODE_DEFAULT_NORM (https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/master/src/common_defaults.jl) both the scalar version and the array version are defined, which is not the case for the provided function u -> norm(u[:, 1]). I think internalnorm should only be applied to compute the norm of the resulting residual, such as in https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/initdt.jl#L126 (which is missing a broadcast as well).

By replacing internalnorm with abs in the computation of residuals I'm able to run

integ = init(prob, alg = Tsit5(), abstol = 1e-6, reltol = 1e-6, internalnorm = u -> norm(u[:, 1]))

successfully and to solve the ODE.

An alternative approach would be that every user has to define a scalar version of internalnorm as well.

[Datseris]

I'd vote for fixing it instead of requiring the user to do more stuff from their side.

[devmotion]

I agree, my only concern was that maybe we should support changing the scalar norm used for computing the residuals as well - but I'm not sure whether this would be of any use to anyone.

[ChrisRackauckas]

but I'm not sure whether this would be of any use to anyone.

Sometimes with new number types you want define a new "norm" there for how to get a value out. That's precisely why it's used in some spots instead of scalar abs. Without this, the only way for someone to input compatibility for their number type would be to hardcode it into DiffEqBase which isn't a nice solution. So somehow we need both norms.

[devmotion]

I fully agree. It seems the best way to handle this is to replace internalnorm with two options for the scalar norm and a norm for u. Then even if users only specify a norm for u (as in the example above) nothing should break.

[ChrisRackauckas]

Added it to the DiffEq v6.0 list since that is a breaking change to how our norms work, though it was never intended at first for the norms to be doubled up like this and just a quirk of history :).

[ChrisRackauckas]

This might cause problems for recursion, vector of vector and other types like that. I think it might be best to just update the docs to say that two dispatches are necessary. Needs discussion of course.

[devmotion]

What are the main problems you're worried about? The default definition of the element-wise norm could still be done in a recursive fashion, it seems to me, so it could also work for recursive types like vector of vectors. I'm wondering whether it might actually be preferable to define the element-wise norm such that abs is applied recursively by default instead of the current implementation in https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/master/src/common_defaults.jl#L2-L7 which computes the vector norm for every element of a vector of vectors instead of applying abs to the each of its scalar entries.

[ChrisRackauckas]

I'm wondering whether it might actually be preferable to define the element-wise norm such that abs is applied recursively by default instead of the current implementation in https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/master/src/common_defaults.jl#L2-L7 which computes the vector norm for every element of a vector of vectors instead of applying abs to the each of its scalar entries.

For the vector of vector case. If it doesn't recursively reduce lower levels to scalars then the result isn't a scalar, and if the result isn't a scalar then the error estimate isn't scalar and adaptivity breaks. We still would need a few other things changed for vecvec to work without VectorOfArray wrappers, but that's at least one requirement.

[devmotion]

For the vector of vector case. If it doesn't recursively reduce lower levels to scalars then the result isn't a scalar, and if the result isn't a scalar then the error estimate isn't scalar and adaptivity breaks. We still would need a few other things changed for vecvec to work without VectorOfArray wrappers, but that's at least one requirement.

Of course, the error estimates have to be scalars, but hence the vector norm should be applied there. I was only thinking about the case of the element-wise norm in lines such as https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/initdt.jl#L125 which shouldn't reduce at all I guess, since the computation of the scalar is only done in the next line https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/initdt.jl#L126.

[Datseris]

So does this work now? Which PR fixed it?

[ChrisRackauckas]

The docs updated.

[Datseris]

I assume you mean this:
internalnorm: The norm function internalnorm(u,t) which error estimates are calculated. Required are two dispatches: one dispatch for the state variable and the other on the elements of the state variable (scalar norm). Defaults are package-dependent.

So this means that for the function f I give to internal norm, I must always define two methods, one for Array and one for Number?

[ChrisRackauckas]

yes