Fixed L2 error for complex valued problems

[henhen724]

Currently, L2 error calculations in DiffEqDevTools.jl use
sqrt(recursive_mean(vecvecapply((x) -> float(x) .^ 2, sol - timeseries_analytic)))
the issue is if the solution is complex this returns the sum of the square of the complex differences instead of the sum of the absolute value squared (i.e. the L2 error should never be complex). For example,

wp = WorkPrecisionSet(prob,abstols,reltols,setups,1//2^(10);numruns=5,names=names,maxiters=1e7,error_estimate=:l2,appxsol_setup=Dict(:alg=>RKMilGeneral(;ii_approx=IICommutative())))
plot(wp)

leads to the error:
isless(::ComplexF64, ::ComplexF64)

I corrected this by replacing x .^2 with conj.(x).*x:
errors[:l2] = sqrt(recursive_mean(vecvecapply((x) -> conj.(float(x)) .* float(x), sol - timeseries_analytic)))

I checked using benchmark tools that this has approximately the same runtime as the line replaced.

julia> function f1!(x)
           return float(x) .^ 2
       end
f1! (generic function with 1 method)

julia> function f2!(x)
           return float(abs.(x)) .^ 2
       end
f2! (generic function with 1 method)

julia> function f3!(x)
           return conj.(float(x)).*float(x)
       end
f3! (generic function with 1 method)

For vec = rand(ComplexF64, 1000), I found the following using benchmark tools:

For vec2 = rand(Float64, 1000), I found the following using benchmark tools:

Checklist

• Appropriate tests were added (I think this change is too minor to necessitate new unit test.)
• Any code changes were done in a way that does not break public API
• All documentation related to code changes were updated
• The new code follows the
  contributor guidelines, in particular the SciML Style Guide and
  COLPRAC.
• Any new documentation only uses public API

Additional context

Add any other context about the problem here.

[ChrisRackauckas]

Add a test?

[ChrisRackauckas]

Could this just be abs2?

[henhen724]

Yes, you could fix it with abs2 instead, but it may be slower (that probably doesn't matter that much). Since the initial pull request I've found new issues. The current version of SciMLBase returns errors that have the same number type as the solution vector. So, for complex solutions it returns ComplexF64 type errors (with zero imaginary part) and then the plotter will refused to plot them. It will also crash while trying to sort the error list because complex number have no isless function. So, I added a function to fix the return type from SciMLBase's function calculate_ensemble_errors and an implementation of isless for complex numbers that just casts them to reals.

The current version has multiple test that are still not working on purpose. Right now, if you provide an SDE problem without an analytic solution you cannot get l2 or linf errors (timeseries errors). So, I am trying to fix that next.

I will send another message when this PR is actually ready.

[ChrisRackauckas]

So, I added a function to fix the return type from SciMLBase's function calculate_ensemble_errors and an implementation of isless for complex numbers that just casts them to reals.

I would think you don't need to change isless and instaed you just want to change the SciMLBase error calculation to always return a real float?

[henhen724]

Yeah, but that would require coordinating two PRs at the same time. So, I thought I better not. Do you want me to open that PR? I already have the fix for SciMLBase.

[ChrisRackauckas]

Well the SciMLBase one needs to happen regardless of anything else, since error should always be a scalar quantity. You might want to just use the diffeq default internal norm for that? Effectively just call norm on the number. I'm curious though what else would need to change, since once you have error as a real scalar I would assume everything else would just work

[henhen724]

Ok, here is the SciMLBase PR: SciML/SciMLBase.jl#787

The other changes I am making here are to add support for non-diagonal noise in some places that it was ignored and to add support for timeseries errors for problems without analytic solutions. Obviously, those are different issues than the one I originally mentioned, but I'm currently using my fork for benchmarking on my main project, so everything sort of got mixed together.

The print statements are from debugging and will be removed in the next commit.