oliverhaas / openabel Goto Github PK
View Code? Open in Web Editor NEWPython implementation of Abel transforms.
License: GNU General Public License v3.0
Python implementation of Abel transforms.
License: GNU General Public License v3.0
Mainly to sort my own thoughts, I will give some remarks on the Hansen-Law method here.
There are two things limiting the accuracy of the Hansen-Law method:
hk
and lamk
.model_hansenLawOrg.hk = [1., 0.596903, 1.09956, 2.57611, 5.65487, 12.2522, 26.0752, 61.5752, 151.739]
model_hansenLawOrg.lamk = [0.0, -2.1, -6.2, -22.4, -92.5, -414.5, -1889.4, -8990.9, -47391.1]
sn[1] = (dataInOld-dataIn[ii])/plan.stepSize
sn[0] = dataIn[ii] - plan.grid[ii]*sn[1]
Typically as of now the method is more limited by the accuracy of the fit. As a note: The interpolation could be fairly easily extended, but would require evaluation of more complicated elementary functions with every order and the interpolation in that order. Both is still linear O(N) computational complexity though.
So a better fit would be nice. Problem is though that the function to fit by exponentials has a singularity at 0 (from here on out Mathematica-like pseudo code):
fun = 1/Sqrt[1 - Exp[-2*x]]
sing = 1/Sqrt[2*x]
As of now I didn't manage to get a good fit. A very robust method to get a good fit is to use some tricks with an SVD to get the exponents first and then do a linear least squares to get the factors. This ignores the singularity, but the whole Hansen-Law method does this anyway. The hope is to get at least a little bit better fit. In Mathematica this looks like this:
nn = 1001;
xmax = 10;
fun = 1/Sqrt[1 - Exp[-2*xmax*x]] - 1;
data = Table[(fun /. x -> z + y), {z, 1./(nn - 1), 1. + 1./(nn - 1),
1./(nn - 1)}, {y, 1./(nn - 1), 1. + 1./(nn - 1), 1./(nn - 1)}] //
Quiet;
svd = SingularValueDecomposition[data];
(*ListLogPlot[svd[[2]]/Max[svd[[2]]]//Diagonal,PlotRange\[Rule]All] \
Enable to manually find sigcut *)
NN = (nn - 1)/2;
sigcut = 16;
lams = 2*NN*
Log[PseudoInverse[svd[[1, 2 ;;, 1 ;; sigcut]]].svd[[1, ;; -2,
1 ;; sigcut]] // Eigenvalues] // Re;
data2 = Table[{x, fun}, {x, 1./(nn - 1), 1. + 1./(nn - 1),
1./(nn - 1)}] // Quiet;
fit = NonlinearModelFit[data2,
Table[a[ii]*Exp[-lams[[ii]]*x], {ii, 1, sigcut}] // Total,
Table[a[ii], {ii, 1, sigcut}], x,
Weights ->
1/Table[fun + 1, {x, 1./(nn - 1), 1. + 1./(nn - 1),
1./(nn - 1)}]^2];
lams/xmax
Table[a[ii], {ii, 1, sigcut}] /. fit["BestFitParameters"]
{326.962, 215.002, 150.651, 108.24, 78.732, 57.6702, 42.4463, 31.378, \
23.3208, 17.4742, 13.2704, 10.2773, 8.02312, 6.00046, 4., 2.}
{4.68558, 1.70955, 1.93147, 1.17423, 1.21348, 0.863059, 0.831268, \
0.635035, 0.578285, 0.453191, 0.386552, 0.304036, 0.284079, 0.312645, \
0.375015, 0.5}
Problem is that we can't chose the number of points nn large enough to get very close to the singularity, since the SVD gets slow really quick. As far as I understand if we're strict, one needs the fit to be accurate in the interval [1/N, Infinity], where N is the number of points in the Abel transform... Too bad overal, because this methods gives very high quality fits/Approximations. In theory one could manually combine this methods for smaller and smaller intervals close to the singularity, but I couldn't manage to make it work in a consistent way.
In principle one could use the standard fitting routines, here I tried using the previous results as starting values among other things:
laminit = lams
ainit = Table[a[ii], {ii, 1, sigcut}] /. fit["BestFitParameters"];
data3 = Table[{Exp[1. s], fun /. x -> Exp[1. s]}, {s, -5, 0,
3./(300 - 1)}] // Quiet
fit2 = NonlinearModelFit[data3,
Table[a[ii]*Exp[-lam[ii]*x], {ii, 1, sigcut}] // Total,
Join[{Table[a[ii], {ii, 1, sigcut}], ainit} //
Transpose, {Table[lam[ii], {ii, 1, sigcut}], laminit} //
Transpose], x]
This converges very poorly. Tinkering around with parameters yields somewhat better fits than the original Hansen-Law, but it's not like there is any consistency in obtaining them, or any real improvement, or any easy way to decide if a small error improvement warrants the use of twice the number of exponentials. And of course the singularity is still a problem, so accuracy is again limit in general.
I tried removing the singularity in the method as well in several different ways. Overall this doesn't seem to really increase accuracy or the ill-posedness of the problem enough.
As of now I'm convinced that the Hansen-Law method is inherently flawed and can't be significantly improved. I'm open for suggestions if anybody reads this :).
For the calculation of the numerical derivative for the backward Abel transform I rely on finite difference. Since the user inputs equidistant data this is fairly natural and works pretty well overall. In general there might be better ways to calculate the derivative numerically (look for Cauchy Integral derivative or at the chebfun Matlab package).
I don't see much value in it now (in my work I can almost always take the derivative analytically), but if there arises a use case one could look into this.
To improve performance there are several places where a "BLAS-based" change of the code could improve performance a little bit. BLAS Cython bindings are available in scipy. As a starting point look at the "convolve" functions used in several modules.
Data interpolation is done with a "polint" function inside several modules, so there are duplicate "polint" functions. Better to move them into one place (probably an interpolate module).
In the same manner the actual interpolation and boundary handling is a duplicate in several modules, better move that to the base module.
Right now the end correction coefficients are loaded and then partly interpolated to the specific transform the user wants to do. The interpolation as of now is cubic, but it would be nice to use higher order for more accuracy if possible and in general not have another duplicate interpolation function (i.e. use "polint" function from point 1)
It might be useful to implement two additional types of transforms. The Abel integral (kernel 1/sqrt(r-y)) and the double Abel transform (which is just applying the Abel transform twice, kernel is just r). I will do this as soon as I need them, which is likely soon, but if I forget I wanted to write down this thought here.
openAbel should be made compatible with Python 3. Since the package mostly is written in Cython, this should be fairly easy to do.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
๐ Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. ๐๐๐
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google โค๏ธ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.