oliverhaas / openabel Goto Github PK

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:

1. Accuracy of the state space model fit of a sum of exponentials to the Abel transform. Hansen and Law have fit a sum of 9 exponentials with accuracy around 1.e-3 when ignoring the singularity. In the code this fit is represented by 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]

2. Order of the assumed interpolation of the input function. Hansen and Law assume linear interpolation of points. In the code this is seen for example by the lines

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.

1. 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).

2. In the same manner the actual interpolation and boundary handling is a duplicate in several modules, better move that to the base module.

3. 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.