Abel transforms of non-noisy & non-experimental data of smooth functions

[oliverhaas]

Hello PyAbel community!

If this is the wrong place to write something like this or if there is no interest at all let me know. I'm still fairly new to the whole collaborating and sharing an github thing.

I stumbled on PyAbel quite a while back and was impressed with the collaboration and documentation on it. Although I couldn't use it directly for my work it helped me quite a bit to get things sorted.

I needed a fast way to calculate accurate forward and backward Abel transforms of large non-experimental & non-noisy data sets of discretized smooth functions and stumbled upon plenty of publications on the topic of Abel transforms. However, either accuracy or run times of the methods weren't sufficient for my work, since the main focus of aforementioned publications is usually the transformation of noisy data and small data sets.

In the end I managed to adapt the Fast Multipole Method (FMM) in combination with high-order end corrections (EC) to my needs and wanted to share some thoughts here. There were similar problems to the Abel transform solved by these methods already published in literature, but since the Abel transform is not one of them up to now and I think there are some pitfalls I managed to avoid I wanted to share.

And if I'm missing something in my explanations or misunderstood something I'm happy to learn.

I implemented the above mentioned FMM+EC and a couple other methods in Python/Cython (repository is here https://github.com/oliverhaas/openAbel, but I just started uploading couple of days ago). Here is an overview of some of the methods, their error, convergence and run times for a test case of a simple Gaussian forward Abel transform:

• HL: Hansen-Law method
• TD: Trapezoidal rule desingularized (same as "direct" in PyAbel I think)
• TE nth: Trapezoidal rule with EC of nth order
• FMM nth: FMM with EC of nth order

One can see that in the case of sufficiently smooth functions the FMM+EC has both linear O(N) run time, high order convergence and small errors. I should note that I intentionally used the forward Abel transform, since it behaves quite a bit nicer than the backward Abel transform. Actually in my work I can circumvent a lot of the problems of the backward Abel transform, because I can differentiate analytically and I don't have experimental noise anyway. But for experimentalists this is of course not the case.

For comparison I tried to do a similar range of plots including both my methods and PyAbel's methods. There are some non-ideal parts of the plots, especially since I haven't dissected PyAbel enough to perfectly call each subroutine, but I think it's a good overview. Legend as above, with PyAbel methods as PA and their PyAbel name. Again a simple Gaussian.

My observation is that accuracy drops quite a bit due to the numerical derivative, even for smooth functions. Convergence order still stays if one takes the correct finite difference of suitable order, obviously. Overall I strongly recommend trying to take the derivative analytically (for the non-experimentalists of course), which in my case was luckily possible due to the convolution derivative theorem and Leibniz's integral rule.

I'm not sure why all PyAbel methods seem to be at most first order. Especially from the Dasch two-point and three-point methods I expected different convergence orders. I'm guessing that most implementations in PyAbel either explicitly or implicitly use first order accurate finite difference derivative.

And I assume the large errors of the Basex method will be mitigated due to the Tikhonov regularization as discussed in #225.

I had some success with FMM+EC and numerical derivatives with maximally flat filters, and if there is interest I can share them at some point. However, typically accuracy is pretty much determined by the amount of noise, and there is often no benefit in going to high order methods. Still, the linear run time of FMM+EC might be useful in some cases as well.

Some special remarks on the Hansen-Law method: I think the approach is quite clever and works well if only low accuracy (~1.e-3) is required. And with that I mean the general approach with the fitting of an exponential model to the Abel transform. I tried to produce fits with more exponentials with only limited success. It's a lot of manual work to get higher accuracy and the end result is slower and less robust than FMM+EC in my experience. Furthermore increasing the order of the integration doesn't really do much either (Hansen-Law assume piecewise linear functions); usually Hansen-Law is limited by the exponential fit accuracy.

Have a nice week.

Cheers

Oliver Haas

[stggh]

@oliverhaas thanks for sharing your nice work!

It would be great if there was a little more information about the FMM+EC method, your particular implementation of Hansen-Law method, and code documentation ;-)

From my quick tests your method (based on your example), it looks good. Here is a transform pair function 3 (1D image size 201):

python2 eg.py
python2 eg.py
openAbel in 4.0 ms, mse=6.56e-10
PyAbel three_point in 5.9 ms, mse=3.54e-10

Is there a pure python version of your code available? Do you have any plans to make available the FFM+EC method for PyAbel?

not sure why all PyAbel methods seem to be at most first order

I think this is because they all employ a similar final step matrix product.

works well if only low accuracy (~1.e-3) is required

PyAbel's most accurate method appears to be three_point

Code for the calculation above. A second run reduces the 3pt execution time to 1ms.

# hacked from example001_gaussian.py (@oliverhaas)
  
import openAbel
import abel
import numpy as np
import matplotlib.pyplot as mpl
import time

# Parameters
nData = 201

# PyAbel transform pair
f = abel.tools.analytical.TransformPair(nData, 3)


# data. This way it's much faster if repeated transforms are done.
t0 = time.time()
abelObj = openAbel.Abel(nData, 1, 0, f.dr)
dataOut = abelObj.execute(f.abel)
t1 = time.time()
mseOA = np.square(dataOut - f.func).mean()
print('openAbel in {:.1f} ms, mse={:.2e}'.format((t1-t0)*1000, mseOA))

t0 = time.time()
pya = abel.dasch.three_point_transform(f.abel, dr=f.dr)
msePA = np.square(pya - f.func).mean()
t1 = time.time()
print('PyAbel three_point in {:.1f} ms, mse={:.2e}'.format((t1-t0)*1000, msePA))

# Plotting
fig, axarr = mpl.subplots()

axarr.plot(f.r, f.func, 'r-', label='analy.')
axarr.plot(f.r, pya, 'g-.', label='PyAbel 3pt mse={:.2g}'.format(msePA))
axarr.plot(f.r, dataOut, 'b:', label='openAbel mse={:.2g}'.format(mseOA))
axarr.set_ylabel('value')
axarr.legend()

fig.suptitle('Inverse Abel Transform {}'.format(f.label), fontsize=16)

mpl.tight_layout()
mpl.subplots_adjust(top=0.87)

mpl.savefig('fmm.png', dpi=75)
mpl.show()

[oliverhaas]

Thank you for the quick response. Many of your questions/remarks could be answered with "Hopefully I get around to it soon to improve it." And I hope that's somewhat normal :). I'll comment anyway.

It would be great if there was a little more information about the FMM+EC method

I'm writing up something more detailed for my PhD, but in principle the references are mentioned in my readme.
The main reference for the implementation of the FMM here was the Chebyshev Interpolation FMM described by Tausch. Tausch solved the Abel integral, which is somewhat related to the Abel transform.
The end corrections are similar to the fairly widely known "Gregory rules" based on the Euler-Maclaurin formula.
These end corrections are described in several publications, but the main reference of openAbel was a paper by Kapur.

your particular implementation of Hansen-Law method

I followed the original Hansen-Law paper fairly closely, but made two smaller adjustments:
For the forward Abel transform I assumed piecewise linear interpolation instead of piecewise constant (or hold) interpolation. And I adjusted the "h_0" coefficients such that it is 1/pi with my definition of pi, since from the derivation that's the one coefficient which should be fixed in my opinion. I didn't have any good reason in hindsight to stray away from the original paper, since the accuracy isn't really getting better in most cases.

code documentation ;-)

If you have anything in particular let me know. I'm improving the code documentation slowly, but if I'm almost exclusively the sole user of it then I want to keep a small (but of course still helpful) amount.

Is there a pure python version of your code available?

Unfortunately not. That's one thing I'm not really planing on doing either. Personally for some calculations Python gets really difficult for me if you need that one single manual for loop and it messes up your whole run time :).

Do you have any plans to make available the FFM+EC method for PyAbel?

Yes I have/had. It depends a little bit on what you're expecting. I planned to write a simple wrapper for sure. But for example if a pure Python version is the minimum requirement for the PyAbel collaborators, then I'd probably stay away from directly contributing to the master branch, and only have my simple wrapper in my own branch and won't bother too much.

I think this is because they all employ a similar final step matrix product.

I'm not sure if I understand correctly what you mean. But I don't think the matrix product has anything to do with it. One could write my FMM+EC in matrix form and it would still be higher order.

PyAbel's most accurate method appears to be three_point
Code for the calculation above. A second run reduces the 3pt execution time to 1ms.

Yes the three_point method seems to have the smallest error for fixed small data sets. Convergence order with increasing N (in PyAbel's case the width of the image in pixels) is still only first order. I actually expected it to have third order convergence since it uses three points and thus a quadratic approximation (EDIT: forgot it takes the derivative, so it should still be second order), but something apparently is not cleanly in the method or implementation, since it only converges in first order.

If I wanted to rewrite your example to be more flattering for my code I would

• Increase the order, which leads to smaller errors, order 5 of your example more or less is machine precision.
• Put the initialization of my abelObj outside the time measurement. Loading coefficients of the end corrections takes a flat amount roughly to do and is done during initialization. Small transforms (less than a thousand points) are not really the best case for FMM+EC, I'm actually happy that it's even close to optimal in those cases. And for my work I do at least several dozen transforms with one initialized abelObj.
• Possibly increase the number of points. Again I know this is not relevant to experimentalists, but if we're talking about transforming analytical functions this is relevant. As I mentioned FMM+EC has O(N) scaling, while all matrix based methods (if the matrices are N x N obviously) have O(N^2) scaling.

openAbel in 0.1 ms, mse=2.06e-13
PyAbel three_point in 0.6 ms, mse=3.54e-10

# hacked from example001_gaussian.py (@oliverhaas)
# @stggh version modfied again by @oliverhaas
  
import openAbel
import abel
import numpy as np
import matplotlib.pyplot as mpl
import time

# Parameters
nData = 201

# PyAbel transform pair
f = abel.tools.analytical.TransformPair(nData, 3)


# data. This way it's much faster if repeated transforms are done.

abelObj = openAbel.Abel(nData, 1, 0, f.dr, order = 3) # default is order = 2
t0 = time.time()
dataOut = abelObj.execute(f.abel)
t1 = time.time()
mseOA = np.square(dataOut - f.func).mean()
print('openAbel in {:.1f} ms, mse={:.2e}'.format((t1-t0)*1000, mseOA))

t0 = time.time()
pya = abel.dasch.three_point_transform(f.abel, dr=f.dr)
t1 = time.time()
msePA = np.square(pya - f.func).mean()
print('PyAbel three_point in {:.1f} ms, mse={:.2e}'.format((t1-t0)*1000, msePA))

# Plotting
fig, axarr = mpl.subplots()

axarr.plot(f.r, f.func, 'r-', label='analy.')
axarr.plot(f.r, pya, 'g-.', label='PyAbel 3pt mse={:.2g}'.format(msePA))
axarr.plot(f.r, dataOut, 'b:', label='openAbel mse={:.2g}'.format(mseOA))
axarr.set_ylabel('value')
axarr.legend()

fig.suptitle('Inverse Abel Transform {}'.format(f.label), fontsize=16)

mpl.tight_layout()
mpl.subplots_adjust(top=0.87)

mpl.savefig('fmm.png', dpi=75)
mpl.show()

2018-09-26 - DH made some small formatting edits.

[DanHickstein]

Hi Oliver!

Thank you so much for showing us your great work on this new method!

But for example if a pure Python version is the minimum requirement for the PyAbel collaborators, then I'd probably stay away from directly contributing to the master branch, and only have my simple wrapper in my own branch and won't bother too much.

PyAbel already has some Cython code (for the Direct method). It would be nicer to keep as much as possible in pure python, but if Cython is required for speed, then that shouldn't keep us from bringing a useful method into PyAbel.

It seems that this method offers access to a very different space of speed and accuracy than our other methods.

The issue of convergence with increasing order is interesting. This is not something that we have investigated much in the past, but I would have thought that the higher order methods would have converged faster as well, and it is strange that they don't seem to.

[oliverhaas]

PyAbel already has some Cython code (for the Direct method). It would be nicer to keep as much as possible in pure python, but if Cython is required for speed, then that shouldn't keep us from bringing a useful method into PyAbel.

A pure Python would probably be nice in general, but I will put my priorities somewhere else in the near future. Of course I would support anyone who wants to give it a shot :).

The issue of convergence with increasing order is interesting. This is not something that we have investigated much in the past, but I would have thought that the higher order methods would have converged faster as well, and it is strange that they don't seem to.

I reread the Dasch paper and as far as I understand it's actually just analytical integration of a piecewise polynomial, so it should be converging with the respective order of the polynomial. I'm guessing there is a small bug somewhere in the implementations.

I actually wrote a small wrapper as mentioned before, so people can use openAbel from PyAbel (see my PyAbel fork https://github.com/oliverhaas/PyAbel and of course you need my code installed as well https://github.com/oliverhaas/openAbel, run/look at example_openAbel.py in the PyAbel examples folder; I'm still working a little bit on both, so don't be surprised). Not a lot of documentation, since it's more or less just for testing/developing.
I think that should be enough tools to decide what and how the PyAbel community wants to progress from here :).

Cheers

Oliver Haas

[oliverhaas]

To add to my previous post, I looked at the comparison of PyAbel and openAbel methods again and tried to figure out the convergence issue. I switched from looking at pointwise convergence to looking at the convergence of the mean squared error (MSE) like @stggh did in his example, which makes more sense in general :). And I changed the test function from a Gaussian to a function which decays to exactly zero at the upper end of the interval.

So here is the result (forgive me for only doing few points for the slower methods):

In combination with my previous results I have some remarks:

• My definition of the "numerical Abel transform" is such that the analytical integral gets truncated, i.e. if R=(N-1)dr then the integral only goes to R and not infinity.
  I think PyAbel might (sometimes?) assume the function decays to zero and appends zero to the input array or sets the last point to zero to calculate the derivative at the last point. This in my opinion produces less consistent results if the user inputs data which doesn't decay to zero.
  In addition above you can see that the errors at the right end of the interval are even large when the function decays to zero, so I'm not sure what exactly PyAbel does there in each method.
• This means two_points works as intended (convergence of MSE and pointwise is 1st order).
• I think onion_peeling probably works as intended as well, but because it implicitly uses a very poor approximation it's only 0th order.
• The method three_points however converges in most points with 1st order, and with 0th order in the MSE. So I'm guessing that there are two bugs in the method: One only at r = 0 and one in general.
• For the basex method it's hard for me to tell what order convergence I should expect, and without the regularization the noise introduced pretty much makes studies impossible. I'm guessing here the convergence depends strongly on the input function, the exact calculations done in the implementation (like truncation of the Gaussians, approximations used, etc.).

A small note: When I say something converges in nth order, it is more correctly (n+1/2)th order due to the square root in the Abel transform.

Cheers

Oliver Haas

EDIT: I replaced the plot with a slightly better one.

[stggh]

1. Sorry to be the dummy here, but what do you mean by convergence and how is it illustrated within your plots?

2. AFAIK the PyAbel Dasch methods make no assumption about the input projection function form. There is no derivative of the input projection image.

3. You are right about the odd behaviour for large R, but for most of the domain the PyAbel Dasch calculations mirror Dash's figure 6:
   

[oliverhaas]

Hello again,

Sorry to be the dummy here, but what do you mean by convergence and how is it illustrated within your plots?

When I talk about convergence in the context of (non-adaptive) numerical methods I mean usually the order of convergence of the relative error of the method with increasing number of grid points for the same problem. For example if you numerically integrate a Gaussian (mu=0, sig=1) from 0 to 1 with trapezoidal rule and use N = 1000 points, you will get a 100 times smaller error than if you would use only N = 100 points. This is because trapezoidal rule is 2nd order accurate (100 = 10^2 times smaller error if you multiply N by 10).
If you plot the error in log-log scale, then you get a just a straight line with the order of convergence as the gradient. You can see in the bottom left plots that this is the case when the method doesn't run into other problems (like machine precision or other sources of errors).
The order of convergence is a good indicator how well an algorithm performs, and if there any bugs in the implementation. In case of the three_point method I'm fairly certain the wrong order of convergence implies some error in the implementation.

AFAIK the PyAbel Dasch methods make no assumption about the input projection function form. There is no derivative of the input projection image.

The method two_point implicitly assumes that the function can be approximated by piecewise linearly connecting the data points, which leads to a constant derivative in between points.
The derivate calculation is a little bit hidden, but it's simple finite difference stencil [1, -1]. Line 154 in https://github.com/PyAbel/PyAbel/blob/master/abel/dasch.py reads

D[Iu, Ju] = J(Iu, Ju) - J(Iu, Ju-1)

Instead one could calculate the derivative of the function first (and make some smaller adjustments) and then just multiply only by J(Iu, Ju).

The method three_point does in principle the same, only with piecewise quadratically connecting the data points. But again, there is an error in the implementation I think. Otherwise it should look different in the convergence plot.

You are right about the odd behaviour for large R, but for most of the domain the PyAbel Dasch calculations mirror Dash's figure 6:

For me it is not clear how Dasch handles the right end of the integration interval. Possibly he doesn't treat it cleanly already in his paper.

Cheers

Oliver Haas

[stggh]

Excellent! Thanks, for the detailed explanation. It had not occurred to me that ΔJ formed the derivative.

[oliverhaas]

Hello again,

so here probably my final remarks:

I looked at three_point again and I think the reason it's only first order is due to Dasch using a first order second derivative stencil at zero. The stencil [1, -2, 1] is only second order accurate if central finite differences are used, which is not done because data is not available at r = -dr. Instead a forward difference stencil has to be used ([2, -5, 4, -1], see https://en.wikipedia.org/wiki/Finite_difference_coefficient) if second order accuracy is desired everywhere. It's probably not relevant for most people, since if only one transform is done this only affects the result at r = 0.

Anyway I rewrote three_point just to check (see https://github.com/oliverhaas/PyAbel) and I get second order convergence (see bottom left plot).

As a note: Please don't use example_dasch_methods.py as a test case, since the combination of how the forward transform is done with hansenlaw and the choice of the data set (non-zero at right end of interval) make it so that by coincidence my version of three_point has a larger error than even onion_peeling.

This will probably the last of my contributions for a while (testing and implementation wise), but I will be around posting and maybe answering questions anyone might have. I know that some of my posts have been a little bit messy :).

Cheers and have a nice weekend.

Oliver Haas

[stggh]

Thanks @oliverhaas, I will take a look at your changes to three_point. You clearly know your stuff!