pbasex

[stggh]

I came across a nice python implementation of the pbasex algorithm:

  https://github.com/e-champenois/CPBASEX

Modifying the provided CPBASEX/examples/sample_pbasex.py provides this O2-ANU1024.txt Abel transform comparison with PyAbel (basex):

test-pbasex.py

There may be a case to work with Elio Champenois to:

1. integrate a version of pbasex into PyAbel
2. improve the algorithm (or parameter option selection) to give a more accurate transform. The rough O2-ANU1024.txt image comparison indicates less resolution, peak intensity, and baseline issues. In particular, the intensity required a xR scaling to match PyAbel, which may even be PyAbel's fault, see abel.tools.vmi.angular_integration() Jacobian should be R^2 not R #253.

Thoughts?

[DanHickstein]

Thanks for finding this @stggh! It would be nice to have a pBasex implementation in PyAbel.

With any new algorithm, we would want to re-make all of the figures from the paper and see how the algorithm compares with the others. In this case, it looks like the resolution is much worse, but perhaps this is simply the basis set.

Please go ahead and contact Elio (if he doesn't already see this thread) and see if he is interested.

[e-champenois]

Hi Dan and Stephen,
Thanks for the message. Feel free to use the pbasex code I developed however you wish, including it in your package if you want. I can of course help out or at least go over what I was thinking when I was implementing this, let me know! I'd be curious to hear more about all these Abel inversion algorithms too, I hadn't considered many of these before diving into pBASEX, since the "pros" made a lot of sense to me at the time.

In terms of the example above, I see two immediate issues:

1. the basis set looks a bit sparse
2. the baseline issue is due to a transpose issue: beta_2 = -1 and beta_2=2 aren't transposes of each other, so the orientation of the image matters with pBASEX.

Here's what I got with your sample image, using 256 radial basis functions (512 total, up to l=2).
I'm sure it could be improved with more tuning of the basis set (number of functions, width of Gaussians, etc.) or other tricks (better centering, l1 regularization, weighting errors, masking weird center noise in the image), but I tried to go with an "out of the box" approach here.

With l2-regularization:

[stggh]

Thanks @e-champenois! Your calculation looks much better. I will attempt to implement pbasex into PyAbel.

[MikhailRyazanov]

I will attempt to implement pbasex into PyAbel.

After #263, please! :–)

... so that I can finish my stuff...

(this is without any regularization)

time in milliseconds

=== Basis generation ==========================================

Method              n = 513    n = 1025    n = 2047    n = 4097
---------------------------------------------------------------
rbasex                114.3       208.7       536.9      2057.1
two_point               3.2        13.9        50.8       192.7


=== Inverse Abel transform ====================================

Method              n = 513    n = 1025    n = 2047    n = 4097
---------------------------------------------------------------
rbasex                  5.4        21.9        83.2       296.3
two_point               1.4        11.3        79.3       671.5

[stggh]

Fair enough!

[haribo3127]

Any update about the implementation of pBasex in PyAbel?

cheers!

[stggh]

My apologies, there is much less free time in retirement than I imagined. I have just returned from sailing in the Australian Whitsunday Islands. I will try to get onto #263 and this task soon!

In the mean time, pbasex is compatible with PyAbel, see my example code above.

[MikhailRyazanov]

@stggh, regarding #263, I'd suggest to work on it two steps: first focusing on the (x, y)/(row, col) issues, and then do the center/origin renaming (which I can do myself, if you don't have time).

[stggh]

which I can do myself, if you don't have time

Sorry to have blocked progress on this for so long. Yes, if you have the time, please start the changes. At the moment I am tied up with documentation, in preparation for a laser safety verification of my lab.

[MikhailRyazanov]

I was reading the pBasex article very carefully, but found that it is not clear about some details, so if anybody can clarify them, it would be great. I'm asking mostly from the perspective of my forthcoming “rBasex” (feel free to try it), but the answers should be also helpful when/if somebody is going to work on implementing pBasex itself.

0. Does anybody have the “official” implementation for testing/comparisons? Are its sources available?

1. Does pBasex really create full images (like in Fig. 2) for each of the M = r_max × l_max basis functions? For our typical images (1 Mpx) this should take about (1000 × 1000 × 8 bytes) × (500 × 2) ≈ 8 GB of memory.

2. In the article they talk about 128 radial and 3 angular terms for 512 × 512 px images. Are these basis parameters adjustable in their program?

3. The article says: “In practice the matrix G is decomposed using SVD into G = USV^T, where U and V are N × M and M × M matrices with N = ij and M = kl, ... The resulting factors are stored on the hard disk”. Since N is the number of pixels (after polar rebinning, but still large, ~million), and M is the number of basis functions (~thousand), the N × M matrix should be also huge (~billion elements, taking gigabytes on the disk and in the RAM). Is this true?

4. After that they say: “The SVD factors can then be read anytime from the hard disk to complete the solution”, and earlier: “The SVD allows us to edit out the singular values to try to improve the condition number of the matrix”. Does this “edit out” mean “remove”/“set to zero” or something more advanced? (The implementation by @e-champenois actually uses an equivalent of Tikhonov regularization, but the wording in the article does not suggest that the original method does.) Are the largest singular values “edited out” (this will definitely improve the condition number), or they are selected somehow else? Can this be adjusted after the singular values are read from the hard disk, or this is done immediately after the SVD step, and the stored values are then used directly in all transforms?

[DanHickstein]

I was reading the pBasex article very carefully

Uh oh. The last time that you read an article very carefully, it was our article, and it was a lot of extra work for the authors! :)

Does anybody have the “official” implementation for testing/comparisons?

I don't have it. I'd suggest that we contact the authors of the pBasex article. I think that their current e-mails are listed here: https://www.synchrotron-soleil.fr/en/beamlines/desirs
and here: https://www.nottingham.ac.uk/chemistry/people/i.powis
I'm happy to have you e-mail them on our behalf, or I can do it if you would prefer.

Does pBasex really create full images (like in Fig. 2)

No idea. Using 8 GB of memory seems fairly crazy...

I also don't know the answers to you other technical questions, but it seems like you are asking important questions.

Do I understand correctly that you rBasex method is similar to the pBasex method, and you're trying to figure out exactly how related the two methods are? Seems like it would be worth getting in touch with the authors of the original paper and figuring out some of the details.

EDIT: the first author of the pBasex paper, Gustavo Garcia, has also done some work on photoelectron distributions from chiral molecules in circularly polarized light: https://www.cfel.de/e4929/e20895/e5101/e5102/2013-12-19GustavoGarcia_eng.pdf
So, if you are contacting him, you may also ask if he has some sample data of asymmetric distributions, as you were requesting in #265 (comment)

[e-champenois]

I can't answer about specifics of the original paper/code but some quick comments:

You are right that you can get into memory issues with large basis sets, though you have an extra factor of 4 in your estimate since there is 4-fold quadrant symmetry in the even-l only case (or at least 2-fold symmetry in the general case, i.e. chiral molecule or w+2w ionization, which could be another nice source of data to test your asymmetric codes: see Zipp et al. Optica 2014).

I ran into these issues trying to fit some above threshold ionization (ATI) data with very high l_max. I think the strength of pBASEX is using the sparsity of the Legendre representation, so perhaps it doesn't make sense to use it anymore (or is not even possible, memory-wise) when this sparsity is lost, which is the case for ATI data.

For practical purposes, N>>M, so the S and V matrices take up very little storage compared to U, which has the same size as G. With U, S, and V calculated, G doesn't need to be saved so the SVD doesn't affect the storage (I guess it might affect the peak memory use unless the SVD can be done in place).

[MikhailRyazanov]

@DanHickstein,

it was a lot of extra work for the authors! :)

Mostly for myself. :–) But in the end everybody's extra work was rewarded, I hope.

Do I understand correctly that you rBasex method is similar to the pBasex method, and you're trying to figure out exactly how related the two methods are?

Yes, since they are closely related, I need to describe what the differences are, but I do not want to make false statements.

I'd suggest that we contact the authors of the pBasex article.

Well, when I asked Hanna Reisler (my PhD adviser) whether we have pBasex, she told me that they did try to get it from the authors, but encountered unreasonable difficulties. Maybe the situations has changed since then, but I still feel somewhat hesitant to contact them. :–) However, if nobody else can help, it's worth trying.

Using 8 GB of memory seems fairly crazy...

I also think so. In the article they talk only about 256×265 (after resampling) images and 128×3 basis functions, which with single-precision floats would give “just” 100 MB, easily doable even in 2004, although still heavy in terms of the disk space. The problem is that this does not scale well to larger images even today...

[MikhailRyazanov]

@e-champenois, thanks for the reply!

though you have an extra factor of 4 in your estimate since there is 4-fold quadrant symmetry in the even-l only case (or at least 2-fold symmetry in the general case

In the original article they do not say anything about the symmetry folding, but even if their implementation actually does it, this is at most a factor of 4, as you said. However, the memory requirements are still in the gigabytes range. And they scale cubically with the image radius.

data with very high l_max.

The scaling with l_max is just linear. I guess, it would be much easier to run out of memory by trying to get a high-resolution transform of a 2048×2048 image even with l_max = 2. So it's not really a sparsity issue.

For practical purposes, N>>M, so the S and V matrices take up very little storage compared to U, which has the same size as G. With U, S, and V calculated, G doesn't need to be saved so the SVD doesn't affect the storage (I guess it might affect the peak memory use unless the SVD can be done in place).

But overall you do need at least one N × M matrix on the disk and in the RAM, correct?

By the way, did you write your implementation just based on the general ideas from the pBasex article, or you had something from the authors? Did you compare with the original pBasex?

[DanHickstein]

in the end everybody's extra work was rewarded, I hope.

Oh, indeed! I'm really happy with the way that manuscript turned out.

I still feel somewhat hesitant to contact them.

Shall I send them an e-mail?

The problem is that this does not scale well to larger images even today

Well, perhaps that's just the way that it works out. Perhaps the assumption is that large basis sets won't need to be used to describe most images...

[MikhailRyazanov]

I actually wrote to Gustavo on Monday and got the source code from him. It is supposed to be compiled by the users and works only in Unix-like systems (he said that other people have created other versions for Windows). I compiled it in Linux without any issues and ran some tests. Yesterday I wrote back to Gustavo to thank him and ask some remaining questions, so I'll report my findings here after he replies. But unfortunately my suspicion about the disk/memory requirements turned out to be correct, and the slowness of the basis generation actually exceeded my expectations...

[MikhailRyazanov]

Now Gustavo wrote that currently he actually uses another implementation, coded in Igor Pro, “with a theta sampling which depends on the radius, which is more efficient and logical”. But since I do not have it yet, here is my summary about the original C implementation.

First, the answers to my questions above:

1. The projected basis functions are indeed created in the RAM. They are not saved to the disk, but instead the SVD is done in place, and the results are then saved. (See below for the memory consumption). The reconstructed-image basis functions (radial Gaussians × angular Legendre polynomials) are computed on the fly (withing about ±5σ from their peaks), so they do not take any extra memory, but require exp() evaluation several times per pixel.

2. The basis parameters are adjustable, but only through editing the source code and recompilation (and then making sure that the saved files match the executable). The number of angular terms is adjustable at run time, but every l_max requires its own basis set.

3. Yes, the matrices will be huge. For some reason the program uses float (32-bit) numbers for some calculations, but the matrices have double (64-bit) elements, so for the default settings with a 256 × 256 resampled image, 128 radial and 2 angular basis functions the U matrix has the size 256 × 256 × 128 × 2 × 8 bytes = 128 MiB (“binary megabytes”) ≈ 134 MB (“decimal megabytes”). Examples for some other sizes are shown below.

4. The program does not modify the singular values at all. The basis-creation part gets them from a library SVD function, which produces M = r_max × l_max singular values (all non-zero, since G is a full-rank matrix), and saves them to one of the basis files. The image-transform part then reads them from that file and passes directly to the library SVD solver. Gustavo wrote: “If you want, you can remove other singular values under a given threshold by setting them to zero”, but the original pBasex implementation definitely does not, and it's not clear whether other implementations do.

Some timing/size tests:

Here NR = radial bins, NTH = angular bins, NK = radial basis functions = NR / 2 (recommended); NL = angular functions = 2 (default) in all cases.

The basis generation has two time-consuming steps: the Abel integration of the basis functions (“basis”) and the singular-value decomposition of the resulting matrix (“SVD”). The decomposition results are saved to separate files for U, V and S; their sizes are added in the table.

NR (NK)	NTH
	128	256	512
128 (64)	6 min basis, 8 s SVD. 16 MB + 128 KB + 1 KB	12 min basis, 17 s SVD. 32 MB + 128 KB + 1 KB
256 (128)		31 min basis, 2.5 min SVD. 128 MB + 512 KB + 2 KB
512 (256)		88 min basis, 20 min SVD. 512 MB + 2 MB + 4 KB	3 h basis, 39 min SVD. 1 GB + 2 MB + 4 KB

So the basis size and generation time scale linearly with NR × NTH × NK (3rd power of the image radius), and the SVD time scales as NR × NTH × NK² (4th power!).

I did not try to create a narrow basis set with NK = NR = 512, which is required to analyze 1-megapixel images with 1-pixel resolution, but it apparently will take >9 h and >2 GB.

The transform timings for a 1 Mpx image are (I have slightly modified the code to print these metrics; the last “Done in ...” is printed by the original code):

Symmetrization: 0.00082 s
To polar: 0.016849 s
SVD: 0.011849 s
PES: 8.6e-05 s
Image: 0.211069 s
Done in 0.241049 seconds

for the default settings (256 × 256 resampled image, 128 functions) and

Symmetrization: 0.000784 s
To polar: 0.060347 s
SVD: 0.094387 s
PES: 0.00012 s
Image: 0.208495 s
Done in 0.365381 seconds

for the 512 × 512 resampled image and 256 functions. Not slow, but not very fast either.

Here I should make a note about the performance of this original implementation. First, its code is not optimized at all, except that it does the left–right symmetrization of the images (but never does the top–bottom symmetrization, even when odd orders are not used) and works with only one half. Otherwise everything is done quite inefficiently (for example, each Legendre polynomial is evaluated at each integration point anew through a recurrent formula, with memory allocation/deallocation...). Second, for all nontrivial calculations it uses the GNU Scientific Library (including the use of its general-purpose integration function to compute the basis projections). The GSL is probably not the most efficient and runs in a single thread (although all the heavy computations can be efficiently parallelized).

I guess that optimizing the computations and using more efficient linear-algebra libraries could improve the timings by a factor of 10–100. The memory usage can be reduced by a factor of 2 (using the top–bottom symmetry) at most. But the time and memory scalings with the image size are still worse than for all other methods, meaning that even with the best optimizations the pBasex method cannot be practically used for high-resolution transforms.

Finally, some additional remarks about this original pBasex implementation and the pBasex article itself.

1. The Cartesian-to-polar resampling is done such that when the image radius is large enough, the NTH angular samples actually skip some image pixels completely. This means that not all the available data is used, and the signal-to-noise ratio at large radii is worse than it could be with proper resampling. Strangely enough, the documentation does not recommend increasing NTH. (I also have an impression that the “bicubic interpolation” is implemented incorrectly.) An implementation with proper rebinning should work better.

2. The comparison with the BASEX results for odd angular orders, where pBasex is claimed to perform significantly better than BASEX, is not really correct. The problem is that the original BASEX program replaces all negative intensities in the output image with zeros (a poor decision!), so when the pBasex authors were extracting radial distributions from that output image, they were integrating all the positive noise (present in the image) without integrating the counterweighting negative noise (removed from the image), and the results were of course badly biased. The BASEX program itself calculates radial distributions without this zeroing. The original pBasex program also does correct calculations internally, but has the same flaw that it zeroes all negative intensities in the output. Thus while the pBasex method indeed produces cleaner reconstructed images, I do not think that its radial distributions are actually better than those obtained by other linear transform methods.

[MikhailRyazanov]

@DanHickstein,

Perhaps the assumption is that large basis sets won't need to be used to describe most images...

It depends on the expectations for what “most” means. :–)

I have shown earlier an example of the 0.3% kinetic-energy resolution (FWHM) in my setup (which was not the first to achieve this), which means the radial resolution ~0.15% and hence at least ~700 distinct radial values, so in practice >1000 basis functions are needed to represent this adequately.

But on the other hand, in many experiments the resolution is indeed much worse and cannot be improved by reasonably better instruments. In these cases a low-resolution method will be enough.

[DanHickstein]

@MikhailRyazanov: That's great that you were able to make contact with Gustavo and awesome that he was kind enough to provide the original source code.

Your analysis of the pBASEX code seem quite insightful!

It seems that, generally speaking, the scaling of the computation time with the size of the image is not especially favorable. How does your rBasex method provide similar results with better scaling?

Your final two points (about the cartesian-to-polar conversion and setting negative values to zero) are perfect examples of how confusion can arise when trying to compare Abel transform methods implemented in different software!

So, have you reached a conclusion regarding if/how we should incorporate pBasex? Or is more study needed?

[MikhailRyazanov]

How does your rBasex method provide similar results with better scaling?

It employs the separability into radial and angular parts (which does not work with Legendre polynomials, but works with cosine powers). So instead of solving the huge problem that connects every pixel to every basis function (where these huge matrices appear), it is splits into two subproblems: 1) get the radial distributions for each angular part (a standard procedure, linear in the number of pixels = quadratic in the radius), 2) transform each of these radial distributions (they are 1-dimensional, so it is a matrix-by-vector multiplication, also quadratic in the radius). The 1st subproblem is as overdetermined as the “huge problem”, but has an explicit solution, without the need of the SVD and extra memory. And the 2nd is not overdetermined, so it's just r_max equations for r_max variables (separately for each angular part).

Plus, the basis projections are calculated analytically, which is much faster, and not for each pixel, but just for each radius.

So, have you reached a conclusion regarding if/how we should incorporate pBasex? Or is more study needed?

It can be implemented directly as described in the article (like the original code does). The only tricky part would be to calculate the basis projections efficiently, if we don't want this to take many hours.

Or, as Gustavo wrote, use a different resampling method (I haven't received his new code yet, but my guess is that it's similar to the rasampling in the “polar onion peeling” method).

And, in my opinion, the approach used by @e-champenois (without resampling) is more reasonable. The only problem is that it is intended to transform images of some particular size, so different image sizes require different basis sets to be generated (which is very slow). The original pBasex resampling just takes its samples from any image, and the basis-set size only determines the relative radial resolution, but not the input/output image size.

[talbert48]

Not certain if this should be made as a new issue or not, but it is sort of in the vein of this thread.

I have been attempting to use the rbasex implementation to invert some above threshold ionization (ATI) data that requires a very high order. I have found that when the order exceeds 20 (or 21 for odds) their is a sharp change in the noise.

The figure below shows an example of this problem. You can produce the same result with the rbasex "beam block" example. The cutoff is always at order 20 (or 21 for odd). As the order increase the problem only gets worse. If this is just a limit of rBasex then it should probably be mentioned in the documentation that their is a maximum order that is usable.

[MikhailRyazanov]

Can you try other “polar” methods (lin-BASEX, pBasex, POP) to see how they behave with very high orders?

And, before going into computational details, is ATI reasonably described in terms of Legendre polynomials? I have no experience with it, but from the few examples that I've read, it seems to be non-perturbative and involve strong interference between many channels.

[DanHickstein]

@talbert48 Thanks for bringing this to our attention, and for posting the nice data.

I think that we should move this to a separate issue. I opened #286 for the discussion.