Inverse Laplace may not be very effectively implemented

[akabla]

I had a look at the inverse Laplace methods we used to see if there is scope for improvement and better integration with expressions. I believe that the Weeks algorithm may be better suited to work on arrays:
https://jlapeyre.github.io/InverseLaplace.jl/latest/index.html#InverseLaplace.Weeks

Here is a bit of benchmark code to compare different approaches (Talbot, Talbot for arrays, Weeks, Weeks in function wrappers, and the current implementation in Rheos based on Talbot in function wrappers). We could go 5-10 times faster than we do now on arrays with 1000 elements, having already increased the precision of the Weeks algorithm compared to the default. As soon as we have 50+ elements, Weeks takes the lead. Weeks is also very convenient to integrate with the function wrappers when an expression is available. We only need to supply the expression of the Laplace transform to the constructor. The constructor would detect that the variable used is :s and not :t and package it accordingly, as done below. We could continue to use Talbot for single scalar values with parameters as it is the most effective one in this case (6 times faster than Weeks for a single evaluation).

using RHEOS
using BenchmarkTools
using InverseLaplace
using FunctionWrappers: FunctionWrapper

t=RheoFloat.(0.0:0.02:1.0);
t0=1.0
cₐ = 1.0; a=0.5; kᵦ = 2.0; kᵧ=3.0
println("Talbot")
Ĵ(s) = (1/s)*(cₐ*s^a + kᵦ)/(cₐ*s^a*kᵦ + kᵧ*(cₐ*s^a + kᵦ))   #FractSLS_Zener
@btime InverseLaplace.talbot(s->(1/s)*(cₐ*s^a + kᵦ)/(cₐ*s^a*kᵦ + kᵧ*(cₐ*s^a + kᵦ)), t0)
@btime InverseLaplace.talbot.(s->(1/s)*(cₐ*s^a + kᵦ)/(cₐ*s^a*kᵦ + kᵧ*(cₐ*s^a + kᵦ)), t)
@btime InverseLaplace.talbotarr(s -> Ĵ(s), t)
println()

println("Weeks")
@btime InverseLaplace.Weeks(Ĵ,80)  # Time to initialize and get the function
ft = InverseLaplace.Weeks(Ĵ,80)
@btime ft(t0)   # single value
@btime (InverseLaplace.Weeks(s->(1/s)*(cₐ*s^a + kᵦ)/(cₐ*s^a*kᵦ + kᵧ*(cₐ*s^a + kᵦ)),80)).(t)

fwp = ((t,p)->(InverseLaplace.Weeks(s->(1/s)*(p[1]*s^p[2] + p[3])/(p[1]*s^p[2]*p[3] + p[4]*(p[1]*s^p[2] + p[3])),80))(t)) |> FunctionWrapper{RheoFloat,Tuple{RheoFloat,Vector{RheoFloat}}}
fwpa = ((t,p)->(InverseLaplace.Weeks(s->(1/s)*(p[1]*s^p[2] + p[3])/(p[1]*s^p[2]*p[3] + p[4]*(p[1]*s^p[2] + p[3])),80)).(t)) |> FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat},Vector{RheoFloat}}}
fw = InverseLaplace.Weeks(s->(1/s)*(1.0*s^0.5 + 2.0)/(1.0*s^0.5*2.0 + 3.0*(1.0*s^0.5 + 2.0)),80) |> FunctionWrapper{RheoFloat,Tuple{RheoFloat}}
fwa = t->fw.(t) |> FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat}}}
println()
println("function wrappers using Weeks")
@btime fwp(t0,[1.0,0.5,2.0,3.0])
@btime fwpa(t,[1.0,0.5,2.0,3.0])
@btime fw(t0)
@btime fwa(t)

println()
println("From the current RheoModel(Class) - function wrappers on Talbot")
@btime creepcomp(FractSLS_Zener, t0, [1.0,0.5,2.0,3.0])
@btime creepcomp(FractSLS_Zener, t, [1.0,0.5,2.0,3.0])
m = RheoModel(FractSLS_Zener,  (cₐ = 1.0, a=0.5, kᵦ = 2.0, kᵧ=3.0))
@btime creepcomp(m, t0)
@btime creepcomp(m, t)

Output:

Talbot
  21.275 μs (642 allocations: 18.75 KiB)
  21.987 ms (642143 allocations: 18.33 MiB)
  3.567 ms (173413 allocations: 3.61 MiB)

Weeks
  116.179 μs (2634 allocations: 110.75 KiB)
  1.610 μs (1 allocation: 16 bytes)
  1.705 ms (2637 allocations: 118.81 KiB)

function wrappers using Weeks
  64.297 μs (71 allocations: 30.80 KiB)
  1.658 ms (73 allocations: 38.84 KiB)
  1.688 μs (1 allocation: 16 bytes)
  1.602 ms (5 allocations: 8.12 KiB)

From the current RheoModel(Class) - function wrappers on Talbot
  10.757 μs (101 allocations: 3.67 KiB)
  10.926 ms (99102 allocations: 3.48 MiB)
  9.649 μs (1 allocation: 16 bytes)
  9.704 ms (1 allocation: 8.00 KiB)

[akabla]

ps: one of the things to figure out is the precision of the Weeks algorithm for our application. It might be that Talbot remains a better option if the numbers are better.

[moustachio-belvedere]

If Weeks is faster, without loss of accuracy, then it's fine with me. But accuracy will need to be evaluated carefully. There is another issue related to InverseLaplace.jl also, maybe there is a way to solve both at the same time. Note I probably don't agree with the solution I proposed there some time ago, I haven't thought about it recently enough to have an informed opinion:

#32

There are benefits to using BigFloat in some of the package's algorithms but I could not figure out how to include that in Rheos in a performant way at that time, and did not prioritise working on it. I'm sure there is something we can do if we want to boost accuracy there:

JuliaMath/InverseLaplace.jl#12

JuliaMath/InverseLaplace.jl#5

[akabla]

Thanks for linking to #32
Do you have a good approach to test accuracy for representative functions?
Each method seems to have its strengths and weaknesses, with an accuracy that may depend on t. Accuracy can also be adjusted in each method by specifying an extra parameter. I think that we should set a default value as a const and allow flexibility to tune for more speed or accuracy.

I think using BigFloat shouldn't be a problem for most things, except NLopt that uses Float64. That is quite a limitation at the moment.

[moustachio-belvedere]

Each method seems to have its strengths and weaknesses

Agreed, definitely worth exploring to find an optimum solution here! And I can see the benefits that weeks brings WRT functionwrappers.

Do you have a good approach to test accuracy for representative functions?

My only thought would be to test the various moduli for which we know both laplace and exact time-domain forms, with some representative parameters, dt and number of samples. But I have not seen any nice documents on truly rigorous testing of complex numerical algorithms.