mpfit.py

[kdere]

It was not necessary to recode mpfit for Python. The same algorithm is available in scipy.optimize.leastsq

hope this is helpful

Ken

[MJWeberg]

Hello Ken,
Good point! Unfortunately, while scipy.optimize.leastsq (or the newer least_squares function) is based on the same minimization algorithm, it does not have all of the same features as MPFIT. In particular, none of the scipy fitting routines allow you to tie the value of one parameter to another, a powerful feature that we use a fair bit with the multi-component fit templates. We also looked at some of the other Python fitting packages, such as lmfit and astropy.modeling, but they either have some operational quirks or are too computationally inefficient for our use at this time.

Thankfully, we did not have to redevelop anything. Although it is not popular within the Python community, mpfit.py has been around and updated since 2002. Additionally, we ran an array of test cases and confirmed that the Python version gives comparable results to the most up-to-date IDL version. Nevertheless, we will monitor the development of other fitting packages and plan to eventually add options for selecting the routine used by eispac.fit_spectra().

Respectfully,
Micah W.

[Richardjmorton]

Hi, do you really need tied parameters when fitting? Scipy optimise minimize offers a couple of options for bounded minimisation (although I do not think tying is available). Plus, you can provide a Jacobian that is calculated using automatic differentiation and combine with JIT compilation (both of which can lead to a speed up of ~x100 compared to the standard minimize speed - I used JAX for something similar). It might be useful considering the amount of spectra that need fitting.

[hpwarren]

Yes, there are a number of cases where tied parameters are used. For example, the centroid of Fe XII 195.179 is assumed to have a fixed offset from Fe XII 195.119. The widths of the Fe XI 188 doublet are another example. They are assumed to be the same, not just within a range.

It seems to me that a spectral fitting package for solar physics needs to be developed. Things that are developed for astronomy are often very slow and don't scale to cases were you have millions of profiles per observation.

[Richardjmorton]

Hmm, from those two examples I think it would be straightforward to avoid using tied parameters (if you wanted to move away from mpfit).

I was wondering a couple of things. a) how long it takes for you to fit a 256 by 256 raster scan with you current method? b) what you time-frame you consider as slow?

While trying out the JAX package I could fit a 256 by 256 scan in ~45 minutes. I am a bit of a novice with JAX, so there could be other ways to increase speed. Plus JAX has GPU/TPU support.

[MJWeberg]

Hello Richard! Sorry for taking so long to get back to you. I am by no means an expert on spectroscopy, but it is my understanding that the tying of certain line parameters is a consequence of the physics (such as line doublets), rather than just a trick for faster fitting.

Thanks for the reference to JAX! I have not experimented with JIT compilers much, but I have heard good things. Currently, we have gotten good results using the built-in multiprocessing library. For example, on a machine with eight cores and 16 GB ram, fitting a double Gaussian to the spectral window containing the Fe XII 195.119 line took ~2.3 min for a 256x256 pixel cutout with tied parameters and ~4.3 min without tied parameters. While the raster I used did contain some off-limb regions (it was the best I could find with those dimensions), this seems like a reasonable range of times based on tests with more complicated regions. It is also worth noting that the data window only has 24 pixels along the spectral axis and that we are scaling the guess parameters before fitting (which gives a moderate speed improvement by itself).

As far as what time-frame we consider “slow”, I am not so sure. I believe there are plans to process large quantities of data eventually, so anything that would make that operation prohibitively long is less than ideal. Nevertheless, quality and continuity of results are more important that just speed and we are always open to ideas for how we can improve the functionality and usability of the code.

[Richardjmorton]

Hi,

If you are interested in this, I been continuing playing with JAX and found it is possible to get much better performance. I managed to reduce the computational time from 45 minutes to 6 minutes - for the same data set mentioned above (256x256x24) on a single Gaussian fit. This was on a CPU with two cores and 12Gb ram.

This would make it possible for nearly anybody to do rapid profile fitting.

Happy to tidy up the notebook and share with you.