Implement fast calculation using hypergeometric function algorithms?

[gonzalezma]

Sylvain Prevost (ILL) would like to implement the fast calculation of 1D and 2D patterns using hypergeometric function algorithms proposed by Wagener and Förster (https://doi.org/10.1038/s41598-023-27558-8), which should be several orders of magnitude faster than standard numerical integration approaches.

Is this feasible within the current architecture? If so, how and which level of effort would be required?

[RichardHeenan]

Myself, Paul Kienzle and Paul Butler had an email discussion which mentioned this towards the end of Feb 2016 whilst investigating the integrations done for exponential shells and linear ramps in the Sasview onion model. However none of this got into a ticket here. I had also at one time started to enquire whether Stephan Forster's code as used in his "scatter" program could be acquired and converted to C++, which would probably be the easiest route. I suspect that to do this from scratch would take quite a lot of effort. See https://www.esrf.fr/UsersAndScience/Experiments/CRG/BM26/SaxsWaxs/DataAnalysis/Scatter

Forster_Burger_hypergeometric_core_shell_eqn_macro1998.pdf

[lucas-wilkins]

I've only looked at this briefly, but I wonder if the comparison to numerical integration is the right one.

As hypergeometric models generalise things like elliptic integrals, Bessel functions, sines and exponentials, any model that uses these is already effectively a hypergeometric model, and under the hood, they will be implemented along the lines described in the paper.

This isn't to say that there are not cases where a more general hypergeometric model would be useful, I'm just not so sure that they will result in much of a speed up in the cases we already cover.

[gonzalezma]

@lucas-wilkins I don't have the answer, but in the paper they claim that they obtain a speed gain of up to 10^5 in the best case. @RichardHeenan I agree that it is likely to represent a significant effort. I have not yet read the paper and even less the accompanying SI (242 pages!), but the interesting thing is that they have already done the conversion to C++ (https://github.com/neutron-simlab/CrystalScatter). The SI also includes all the equations and the code in Mathematica.

[lucas-wilkins]

The SI also includes all the equations and the code in Mathematica.

Hehe, as soon as I saw the words "Hypergeometric function" I thought "sounds like someone's being doing integrals in Mathematica"

[butlerpd]

LOL .. so as it turns out, I have been trying to get my hands on SCATTER for several months now. However, Stefan Forster was on our panel at NIST recently and so I had occasion to discuss this with him directly. I was pleased to find out that he'd found somebody to translate his old Pascal code into C++ and he indicated they were working on getting a real release of this new version in the next few months? I've been watching the github site but not much activity there so far.

Stefan was certainly supportive of having SasView implement some of these and apparently he'd already been approached by Toby Perring at ISIS about whether that would be feasible?

From my reading of the paper and discussions with Stefan, my take was that we should be able to gain some speedup for the models he's done the math for (appendix), though It may just be the numerical integrals from what @lucas-wilkins says (not sure I understand at that level TBH). However I don't think we can get the full speedup he gets in any case. The place we can "easily" start is by rewriting (or adding an alternate method) for form factors he's converted to hypergeometric functions. In fact that could be a good place to start the benchmarking of what speedups we can get?

On the other hand, polydispersity integral won't work at this time (and not sure will ever work for more than specific types of polydispersity - he uses the Shulz-Zim which is great for colloids) and the structure factor integrations he does I also don't think are easily translatable? He also gets speedup in this latest work (which as I understand it was not in the original SCATTER) by pulling out as much as he can into a q independent calculation which then gets calculated once rather than each q. All of these would require changes to sasmodels and therefore some careful design work? Also he currently does not support non-centro symmetric particles like Janus particles which require complex math. It is not clear to me that would be possible but my impression was that he thinks it might be?

Still, as I say, I've been thinking it might be worth writing the general parallelepiped for example to see how the 2D computation speed improves? Just a thought.

By the way @RichardHeenan, I originally remembered our conversation as you did and looked for a ticket, but then remembered we had a discussion at one of the Grenoble code camps. Turns out it was the last code camp before the pandemic and there was an action item to revisit the following year https://github.com/SasView/sasview/wiki/CodeCampIX_BreakOut1 😄

[lucas-wilkins]

I've been listening in to the ML workshop at Berkeley, someone just asked this exact question (before I did). The answer was its not going to be any faster than e.g. Bessel functions.

[pkienzle]

2D parallelepiped in particular is fast. It is just three sinc functions. Hollow forms require two sets of sincs (interior and shell). Using linearity skew transforms can be applied to Q before computing the sinc. Are there other generalizations that you have in mind?

There may be speedup if you can use analytic forms for angular dispersion, size dispersity and resolution. Otherwise we are doing a gridded implementation of a high dimensional integral, which quickly gets slow as dimensions are added.

More generally we could switch to quasi Monte Carlo integration to speed up all the models (see #132).

[lucas-wilkins]

More generally we could switch to quasi Monte Carlo integration to speed up all the models (see https://github.com/SasView/sasmodels/issues/132).

I was actually going to implement something along these lines next week. Might be good to chat with you if that's something you had in mind.