-
Notifications
You must be signed in to change notification settings - Fork 8
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Bromwhich contour: Hyperbola; comments on inversion at t = 0 and fixed contours #23
base: master
Are you sure you want to change the base?
Conversation
Trying to sync with main package now
Thank you, this seems like a very interesting and valuable PR and I appreciate the highly informative PR message. Thanks also for keeping us in the loop. JuliaRheology/RHEOS.jl#32 is certainly still a topic of interest for us, though work on RHEOS is now somewhat intermittent. With regards to how the functions should look, there is an example in #12 , the time domain form is
I don't have to hand all the other functions we tested in time domain + Laplace domain but most of them rely on either summed exponentials (similar to example above) or some power law terms. So for another test case this might be a good one (hopefully I have both forms right):
Lastly, I would just note that on JuliaRheology/RHEOS.jl#110 I don't believe the title refers to any bad implementation of InverseLaplace.jl itself, rather that us RHEOS developers may have not integrated it into RHEOS.jl in the best way, and the interest in exploring the Weeks algorithm which may be faster in speed-critical part of RHEOS. I hope the above is helpful, and thanks again! |
Sorry I missed the important part:
That sounds right to me. From physical principles I would expect it to certainly go to 0 as t->0, and if it didn't I would assume I had made a mistake in the analytical form. Unfortunately this was one of the moduli we were not able to analytically invert so I don't have more insight. |
Thanks so much for the comments!! And very sorry for the long delay as I got busy writing proposals. These two examples you provided are excellent as they are so different, and maybe this can help with your discussion on JuliaRheology/RHEOS.jl#110. I'm going to try to discuss this in terms of 1. accuracy and 2. performance. AccuracyI attached the salient parts of the code I ran to produce the plots below. ## two function to look at
f1(t) = t^0.5
g(s) = gamma(1.5)/s^1.5
f1(t) = 1 - exp(-t)
g(s) = 1.0 / s - 1.0 / (s + 1.0)
## compute Weeks function weights, performance critical
w = Weeks(s -> g(s), 80)
## define our t array, the ratio of t_end/t_start affects accuracy of fixed hyperbolas
t = range(0.05, 40.0, length = 60)
## compute absolute difference for hyperbola, talbot, fixed hyperbola, and weeks
## hyperbola, talbot, weeks take a float value, hyperbola fixed requires an array
h_err = hyperbola.(s -> g(s), t) .- f1.(t) # hyperbola
t_err = talbot.(s -> g(s), t) .- f1.(t) # talbot
w_err = w.(t) .- f1.(t) # weeks
hf_err46 = hyper_fixed(s -> g(s), t, N = 46) .- f1.(t) # fixed hyperbola using 46 evaluations
hf_err76 = hyper_fixed(s -> g(s), t, N = 76) .- f1.(t) # fixed hyperbola using 76 evaluations For these two examples, the hyperbola contour converges the quickest and is the most accurate. There is a small caveat that the Weeks method produces surprisingly accurate results for Based on these examples and others as well. The hyperbola contour is able to converge quicker and more accurately than the talbot contour. Also interesting is the ability of the fixed hyperbola to match the accuracy of talbot given enough evaluations. In general the larger the ratio of t[end]/t[1] the more evaluations will be needed. I showed here 46 and 76 evaluations. Depending on your application, the 46 evaluations can give about 7 digits of accuracy. Again for a smaller time window, less evaluations would be needed for similar accuracy. It is very strange to note the consistency in the error in the fixed_hyper for both evaluations... I honestly, am not for sure why the error is so stable. It is quite interesting. The advantage of this approach will be highlighted in the performance section. Performancet = range(0.05, 40.0, length = 2000)
@btime hyperbola.(s -> g(s), $t)
9.283 ms (1 allocation: 15.75 KiB)
@btime talbot.(s -> g(s), $t)
10.472 ms (1 allocation: 15.75 KiB)
@btime w.($t)
2.212 ms (3 allocations: 15.88 KiB)
@btime hyper_fixed(s -> g(s), $t, N = 46)
1.866 ms (22 allocations: 18.59 KiB)
@btime hyper_fixed(s -> g(s), $t, N = 76)
3.152 ms (22 allocations: 19.66 KiB) In order of increasing computation time: There is this snag when |
Dear Michael @heltonmc , |
It looks like I somehow missed this PR. It's a year old! Thanks for doing this. I do have a concern about the many floating point literals and conversions in this PR. I'm not certain that they matter. For the Talbot routine they certainly do when you use For example And in the examples, it's better to use |
But, it makes sense that for some applications |
No problem! And no worries if this isn't wanted in this package, I should have filed an issue first. It provides similar features to With that said reviewing this code with a bit more julia experience, I did have a couple concerns I'll mark if you want to go ahead with this PR which I can fix. |
src/bromwhich_contour.jl
Outdated
""" | ||
function hyperbola(f::Function, t::AbstractFloat; N = 16) | ||
a = zero(Complex{eltype(t)}) | ||
N = convert(eltype(t), N) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Probably not needed to convert types here as it will be type stable without this if N
is an Int. We should limit N
to be an integer as line 57 will not be accurate if it's not.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes. For a number of reasons, you don't want to convert N
.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes. For a number of reasons, you don't want to convert N
.
src/bromwhich_contour.jl
Outdated
|
||
# parameterized hyperbola contour | ||
# parameters are given by ref [2] which are given in increased precision from calculations given in [1] | ||
function s(θ, N, t) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
s
and ds
should be combined in a single functions as they share several arithmetic operations. Can save at least one multiplication and division per loop doing this.
I don't think there's an option to return both sinh
and cosh
as there is a sincos(x)
function.
src/bromwhich_contour.jl
Outdated
|
||
# get the fixed integration components | ||
function hyper_coef(N, t::AbstractArray; ϕ = 1.09) | ||
A = acosh(((π - 2 * ϕ) * t[end] / t[1] + 4 * ϕ - π) / ((4 * ϕ - π) * sin(ϕ))) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
O no!!! 😲 Considering all the recent blogs and attention on this t[1]
seems like a potential problem here considering we allow t
to be an AbstractArray
. This also assumes t
is ordered and t[1]
is the start or first value in t
.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes, but it's really easy to do this. I think first(t)
will give you what you want.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
But, I'm pretty sure AbstractArray
has to be indexable, and the indices fix an order.
src/bromwhich_contour.jl
Outdated
dsk = ds((k + 1/2) * h, N, t) | ||
a += f(sk) * exp(sk * t) * dsk | ||
end | ||
return imag(a) * h / pi |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I use pi here and π
below which is a style inconsistency.
src/bromwhich_contour.jl
Outdated
# compute the function values over the fixed contour nodes | ||
function fixed_sk(f::Function, N, t::AbstractArray) | ||
μ, h = hyper_coef(N, t) | ||
a = zeros(Complex{eltype(h)}, Int(N)) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Hmmm I always go back and forth on the best way to handle this. But we are allocating quite a bit within this function that the user can't do anything about it. And then we call Threads
which greatly increases the allocations as well. If the user is calling this within a loop they'll see a lot of allocations. Wondering if she provide some pre cache option here....
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
You should probably used typeof
rather than eltype
for things that aren't containers. Although eltype
works. Maybe because floats and ints are iterable.
It's common to provide fixed_sk!
that takes a buffer as input and then fixed_sk
just does this allocation and calls fixed_sk!
.
I don't know the best way to handle threads; you can provide a keyword argument. You can also have a conditional that turns on threads only if the work is big enough.
If I understand, these routines provide some significant improvement over other methods for some problems, right? Converge more quickly, or with better accuracy? If so, then I think we should go ahead with this. I am traveling for the next few days. But, ping me if you don't hear from me. |
@jlapeyre sorry for the delay. I think I've now fixed all these issues and ready to review. All the functions will now preserve input types so you can use these functions with But yes, these are improved versions of the original |
Hyperbola
This pull requests implements the numerical inversion of the Laplace transform by using a hyperbola for the Bromwhich contour. It's important to say that the inversion is very ill posed so different contours will lead to different convergence and accuracy.
Major Changes:
bromwhich_contour.jl
was added that utilizes a hyperbola for the Bromwhich contour.hyperbola
uses a contour that depends on the time point which yields most accurate results.hyper_fixed
computes a single optimal contour for a given range of times needed that is much more efficient at the cost of accuracy.hyperbola
andhyper_fixed
.Comparison to Talbot contour
As an example, let's compare the convergence of the Talbot and hyperbola contours for
![Screen Shot 2021-05-06 at 8 11 19 PM](https://user-images.githubusercontent.com/57471570/117381403-fcdd4e00-aea9-11eb-9468-e67ce854b154.png)
t = 1.0 and 10.0
for the functionsF(s) = 1/s - 1/(s + 1), f(t) = 1 - exp(-t)
.Take aways:
Inversion at t = 0
I have seen a couple issues regarding this specifically #12 and #5 which has led to some confusion with packages using this JuliaRheology/RHEOS.jl#32 and JuliaRheology/RHEOS.jl#110.
The integrand we are trying to calculate looks like
exp(st)*F(s)ds
. Theexp(st)
term is highly oscillatory on the Bromwhich line andF(s)
typically exhibits a rather slow decay. So now we have an integral that is both highly oscillatory and decays slowly..... this is hard to compute. The ingenious idea by talbot was to circumvent this slow decay by deforming the contour to begin and end in the left plane, i.e.,Re s -> -inf
at each end. Now,exp(st)
when evaluated ats -> -inf
forces a rapid decay of the integrand so we can then approximate the integral by the midpoint or trapezoidal rule.Now, if
t=0
we see that the exponential term no longer forces rapid decaying of the integrand so we must evaluate larger values and more terms.Alternatively:
From the optimization approach we are trying to determine the contour parameters that give fastest convergence. We choose these parameters to be proportional to
N/t
. In fact if you look at the source code you will see:So if
t = 0
you get one of the parameters to be infinity which leads to undefined behavior. Thus, the original Talbot implementation was operating appropriately, it just can't be evaluated whent = 0
. **Note:talbotarr
works whent[1] = 0
because the contour is fixed byr = (2 * bM) / (5 * tmax)
, so we never end up dividing by zero. This approach though is less accurate though I have not extensively used fixed talbot. The original suggestion was to approximate the value att=0
by evaluating at a small positive value close to0
. Let's look at the contour we are integrating in the real and positive imaginary plane as a function of time. (only need to consider half of contour for real valued functions).If you have really good vision you can see the blue curve when
t = 1.0
is a tiny dot near 0. So we can see that the contour size depends on t which is important for the evaluation near0
as well as the topic I'll discuss later considering a fixed path. For values close to zero we have to accommodate the larger contour with more function evaluations.Solution
As originally shown by @jlapeyre, we can't directly evaluate at
t = 0
, but we can approximate att = 0
by evaluating very close to0
. In #5, we had this original example:I'm not for sure what the analytical
f(t)
looks like here, but it appears to go to0
ast -> 0
and towards1
ast -> inf
. Can @moustachio-belvedere comment on how this function should look? But let's quickly look at a case when we know the solution as shown in #12. Again this is expected behavior as discussed previously, and I believe the work around is being implemented correctly and further resolves #5 and resolves #12.Evaluating for many time points
So I wanted to try and discuss different issues originally brought up in #5 but seem to be open in issues JuliaRheology/RHEOS.jl#32 and JuliaRheology/RHEOS.jl#110.
First,
talbot
andhyperbola
requireN
evaluations of the functionF(s)
for each time point. The parameters of the contour depend ont
, so for each new value oft
,F(s)
needs to be recomputed on a different contour. This is quite slow whenF(s)
is computationally expensive or many time points are needed. As I've discussed and shown in the previous figure, the contour size and shape depends ont
. In fact, it is optimized for each individualt
point for best convergence.However, this has a drawback because for each new time point we must also compute a different contour and evaluate our function again. What
talbotarr
andhyper_fixed
seek to do is a pick one contour shape that does not depend ont
(only the ratio oft[end]
andt[1]
and evaluate the function onlyN
times instead ofN*length(t)
times. So now this method is significantly faster proportional to however long it takes to computeF(s)
and how many time points needed. For one of my applications it is over 5 orders of magnitude faster.The question now is how accurate do you need the calculation to be? Instead of having an optimized contour for each time point we now have to compromise to an average best contour which will decrease accuracy. The decrease in accuracy will depend on the ratio between the start and end times. Let us try to show us an example of what I mean:
Ok so we can see that if the range of the start and end t values is too large that increases the size of our contour and will require more evaluations. To resolve JuliaRheology/RHEOS.jl#110, I am not for sure what the true values should be so I can't comment on the best way to invert this. @akabla or @moustachio-belvedere may have to provide expected values. Though using
talbotarr
as well ashyper_fixed
will provide poor results if your start value is very close to zero and your end time value is very far. Going back to the previous example:We can increase the number of evaluations so this test passes even with a really low start time. Honestly, this is pretty incredible that this test past when
t[end]/t[1]
is100000
. This approach will be orders of magnitude faster than usinghyperbola
ortalbot
. I'm hoping some of these explanations will help close JuliaRheology/RHEOS.jl#32 as well. One thing to mention that would be a hybrid approach is to utilize 2 contours. Usehyperbola
to estimate a value close to zero and then usehyper_fixed
for the contour so that the range of time values are smaller. This could actually allow for less computations ifhyperbola
converges very fast then by decreasing thet[end]/t[1]
ratio you would need less evaluations there as well. Though this really depends on how close to zero you need to knowf(t)
. If 0.1 is sufficient probably can just use the same contour but if we need1e-100
.. then might be best to consider two contours.@jlapeyre There are also seems to be some confusion about some of the implementations of the approaches I've seen through the issues. I would have some time to work on documentation more extensively but wanted to know if you had a vision for this, maybe we can open an issue to discuss separately.