Good news for all of you keeping up with current affairs in the worlds of Technology, Entertainment and Design: There is a TED conference coming to Europe this summer! So if you are unwilling to travel to the US, you can now travel to the UK instead, and join the conference from 12-16th of July in Oxford. Registration and information available here: http://conferences.ted.com/TEDGlobal2010/

I probably will not be able to convince my boss that this is a conference I need to attend, and I’ll be in Japan that time, so I think I may unfortunately have to miss it due to a lack of funds. But for all of you that are going, let me know how it was!

With the help of Michael Croucher of the walkingrandomly weblog, who graciously spent some time optimising one of the slowest components of the Matlab code, we were able to come up with an approximation to the Bessel function which was about two times faster than Matlab’s internal function. However, then I had to smoothen out the crossover point using error functions, and now it has lost some of the speed it had gained in favour of smoothness. It is now about 40% faster than Matlab’s Bessel function, when approximated using only 5 polynomial terms.

It looks like this, with the crossover gradient centre at 4.5:

So there remains a small discrepancy for such low polynomial numbers, but it should be possible to obtain a good speed increase with this in cylinder and disk scattering functions. Please grab yourself a copy of the code, and I’ll be happy (albeit perhaps slow in the coming two weeks) to include more improvements!

bessel_li_gross

p.s.: I forgot to mention how to compile the code for your computer. At this point, there’s a compiled version available for intel Macs. If you do not have that one, but you do have the matlab (”mex”) compiler, I suggest you follow the following steps:

1. make sure your mex compiler can find the following include files: mex.h, math.h, omp.h (in my case, I copied a suitable version of omp.h into the bessel_li_gross directory

2. go to the bessel_li_gross directory and type “mex blg_component.c”, or if your compiler complains about not being able to find the omp.h file that you just put in the directory: “mex -I./ blg_component.c” (adds the current path to the search path for files under unix)

3. Play. You can also send me your .mex file so I can add it to the package.

The Bessel function is used in a few scattering functions, for example in that of a cylinder. The implementation of the Bessel function in Matlab is quite slow, however, which annoyed me considerably as it slowed down my fitting functions. I therefore had a look at whether I could speed it up.

The paper by Frank B. Gross entitled “New approximations to J0 and J1 Bessel functions” (IEEE transactions on antennas and propagation vol. 43, no. 8, 1995, pp.904) seemed like a good place to start. It is mentioned there that there are two fast approximations already available, i.e. the “small argument approximation”, valid for x<0.5 and the “large argument approximation”, valid for x>5, which are very quick functions indeed. However, for 0.5<x<5, the only approximation available was a Taylor series.

dr. Gross then comes up with a new approximation expressed as a polynomial with integer coefficients. The only drawback of this implementation, is that for the first order Bessel function J1, it requires the calculation of the derivative at x, and would therefore require either more points, or would put stringent requirements on the order and spacing of x.

In a later paper by Millane and Eads (IEEE transactions on antennas and propagation, 2003, vol. 51, no. 6, pp. 1398) the approximation derived by Gross is further generalised to encompass all Bessel functions of integer order. More interestingly, it does away with the derivation. However, I could not get their approximation to result in the correct values.

The last paper that “solved” the issue is a paper from 2006 by Li (and Gross, Applied mathematics and computation 2006, 183, pp. 1220), which generalises the generalisation by Millane and Eads to obtain an approximation for the Bessel functions for all fractional order. This contains an equation which provides a good and quick approximation to the Bessel function.

When implemented, however, it could not beat the internal function in Matlab. The speed reached by the approximation was virtually identical to that of the internal Matlab function. Now, while this was not a great success, it can be considered a partial success, as my interpreted code was as fast as the optimized, internal code. Of course, it remains an approximation…

Finally, I tried simply using a linear interpolation between several results of the Bessel function between 0<x<5, and using the large argument approximation beyond that. This gives moderately good results, as it is two to three times faster as the internal Matlab function. While not much faster, it has certainly taught me a lot about Bessel function implementations. If there is anyone among you out there with knowledge on how to compile Matlab functions into optimized code, please let me know.

Cheers,

Brian.

During the past few months I have defended my thesis successfully and I have started work as a post-doc at Aarhus University working on SAXS studies of reactions under supercritical conditions. I will update this website in 2010 as well, so keep tuning in!

I wish everyone a very successful (and fun) 2010!

The two binning routines made available bin _only_, they do not normalise to the number of pixels in the bins. In order to achieve this, adapt the weighted binning routine by changing line 39 to:

Ibin(bini)=sum(Isort(qi(1):qi(2)).*weightfact)./sum(weightfact);

I have not looked at a similar solution for the equal intensity binning routine yet.

Hello,

Once again, I’d like to present a binning routine. This time it’s a weighted binning routine, which instead of simply summing all intensity in one bin of a certain width, it splits the intensity to neighbouring bins based on their proximity to the bin centres:

As is explained in the graphics, this implies that there is a loss of intensity required in order to keep the first and last bin in line with the other bins. The advantage of this binning function, however, is that it should result in a smooth binned curve with no moire-like effects.

Speedwise, it bins a 2048×2048 image in about 7 seconds on a 2.4 GHz Intel core 2 duo processor (single core), which is just as quick as the equal intensity binning routine discussed in a previous post. As always, if you have any suggestions or improvements, please let me know so I can update the code.

The code itself can be found here: binning_weighted code

It is also given below (behind the “more” button). Leave comments please!

 Read the rest of this entry »

The famous question of the uninitiated in small-angle scattering is: “Do you have a bit of software which will give me an answer from my data?”. After a lengthy explanation (coloured with some anecdotes) about why small-angle scattering is not a uniquely defined problem with an often unique answer such as wide-angle diffraction might be, the new user is then left with a copy of Matlab or Fit2D and asked not to return until he has a more “sensible” question. I guess this is because there are not many alternative treatments for these users. These days, I may also give them a copy of the most recent SAXSGUI, but since this lacks quite some fitting functions, it is only really useful for users who already know what they are doing and can program their own 1D or 2D fitting functions. An all-in-one package that is not only good for beginning users but can also remain a useful tool for advanced users is not something which I’ve seen so far. Until now, that is.

I’d like to draw your attention to a very useful software package for anyone working with isotropic small-angle x-ray/neutron scattering data called SASFIT and can be found here. I noticed this package only last month during the SAS09 conference in Oxford, but I have been told that it has been under development for quite some years now by Joachim Kohlbrecher and Ingo Bressler. It comes precompiled for OS X (Macintosh) as well as for some other unmentionable operating systems, so colour me happy.

The one thing that struck me as particularly good about this package is its completeness. I think you will be hard pressed to find any other package containing as many form factors and structure factors to fit your 1D data as this package. On top of that, all the functions are described in the documentation, so if nothing else, do take a look at that (and ste– recode the functions in your favourite language of choice if you must).

The code is open source (as opposed to some other SAS fitting software which is more like a black box) and it has been written in C. While it is not the language I’m familiar with, it is a very good choice. Much more readable than FORTRAN (in which much of the other software is written), can be compiled on a plethora of platforms, and is a language which is the backbone of much good software (i.e. it doesn’t change much over time, it has been around for a while, and it will be around for a while).

To get going, you must have your data in a format readable by SASFIT, and the easiest data format is to have your data in an ASCII-readable file with three columns, q, I, and the error in I (although the last is not necessary). That means you must preprocess your data using your own tools to obtain it in this form, for example using the binning method described two weeks ago, which will supply you with data-points with equal error. Once that is done, you can use your binned data as input into SASFIT.

When using the software, I must say it was not the most intuitive software I have ever used. The TCL/TK interface is partly to blame for this, as it’s not very well integrated with OS X, and therefore appears awkward and occasionally “hangs”. Reading the manual is therefore recommended for this software, to prevent premature baldness. Once you get going, though, it’s very good. Initially, applying some random fitting functions to my data resulted in many (not very straightforward) errors, which could be resolved when setting the parameter boundaries for the fitting parameters correctly. If you do not do that the fitting may optimise parameters to “impossible” values, which subsequently results in impossible intensities, impossible residuals, and impossible optimisation criteria. Once you set the limits correctly, it works much better, and good results can be obtained.

During my initial test, I found the optimisation algorithm to be not so robust. This is likely due to my choice of fitting function, but from my experience with the Matlab “fminsearch” and “fminsearchbnd” functions, I had expected these to work a little better. The fitting is also rather slow, which may be due to the continuous updating of the fitting curve in the graphical display.

These experiences of an inexperienced user aside, I would very much recommend this package to anyone working with 1D SAS data. It holds much promise, contains a hell of a lot of fitting functions, is well documented and open source, and in my eyes is essential software for all levels of users. If it becomes slightly easier to use off-the-bat (i.e. with preset sensible parameter limits for each fitting function, some pre-loaded example data and a good welcome screen pointing you at the place to start) our lives would become a lot easier: “Take this package, little man, and don’t come back until you have a sensible answer!”.

Hello everyone,

I’d like to draw your attention to a course I would recommend to any ph.d. student using neutrons or X-rays. I’ve had good fun during this course, and learned a lot, so I’d recommend anyone to apply for this:

Read the rest of this entry »

Binning of data is often done to SAXS data in order to reduce the size of the data, “improve statistics” or otherwise retransform the data to better suit the subsequent plotting or fitting of the data. Common 2D and 1D binning methods are “linear” binning, where the q-range is divided into bins of equal size, and “logarithmic” binning, where the q range is divided just so, that the bins are equidistant on a logarithmic q-scale.

Depending on your binning method, you may be weighting your data differently in subsequent fits of no weighting is done by uncertainty or error of the intensity. If you would like to weigh your data points in the fitting by their errors, you need to compute the error for each pixel, which is often given as the square root of the number of counts on your detector. This does imply that you either need a single photon counting detector or that you need to scale your detector output to photon counts.

An alternative to this is to rebin your data just so, that all points have the same intensity, meaning that you have wider bins at higher q than at low q. Since each bin then contains the same relative amount of intensity, the error in each bin should be the same (similar solutions can be achieved with non-position sensitive detectors such as the cyberstar, by counting in each position until a certain threshold value has been reached).

I have written a small Matlab program which achieves just this, and is fast enough to do this for large datasets (i.e. 2D images consisting of 2048×2048 pixels) within seconds. I invite you all to take a look at it, and, as always, suggestions or alternative methods are more than welcome. The input is self-explanatory, but the output isn’t necessarily. The output consists of four vectors:

• qbin: the mean centres of the bins
• dqbin: the width of the bins
• Ibin: the intensity in each bin
• Ibdq: the binned intensity divided by the width of the bins. This is given for plotting purposes (since plotting qbin vs. Ibin would give you more or less a straight line).

Good luck!

equal intensity binning routine:

 Read the rest of this entry »

After too much time, I have finally upgraded to a newer version of wordpress. Previously I had some issues during writing, so that it did not easily understand paragraph breaks, I hope that’s finished now.

In other news, I am only 9 days away from handing in my thesis. The name of the website may change after that, but since I’ll still be working on SAXS, it’ll certainly continue. If you have any suggestions for content or just to say hi; drop a comment below!

B.