[SciPy-dev] Mpfit optimizer

[Bill Flynn]

And I realize that I forgot to attach our files. We have them in our svn
repo at
svn checkout *http*://
spinwaves.googlecode.com/svn/trunk/spinwaves/utilities/mpfit

We have the mpfit.py file and some tests for it as well.

Thanks!

Bill Flynn


On Wed, Aug 19, 2009 at 4:09 PM, Bill Flynn <wflynny at gmail.com> wrote:

> Hi,
>
> I received an old python version of the mpfit optimizer earlier this
> summer. My advisor had previously been in contact with SciPy about including
> this optimizer in the optimize package. He was sent to following message
> with tips to help improve the package:
>
> """
>
>> ---------- Forwarded message ----------
>
> From: Pauli Virtanen <pav at iki.fi>
>
> Date: Sat, May 9, 2009 at 11:55 AM
>
> Subject: Re: [SciPy-dev] adding mpfit to scipy optimize (Please respond)
>
> To: scipy-dev at scipy.org
>
>
>>
>> Fri, 08 May 2009 16:09:35 -0400, william ratcliff wrote:
>
>
>> > Hi!  For a long time, there has been a standard package used by IDL
>
> > users for fitting functions to data called MPFIT:
>
> > http://cow.physics.wisc.edu/~craigm/idl/fitting.html
>
> > <http://cow.physics.wisc.edu/%7Ecraigm/idl/fitting.html>
>
> [clip]
>
> > http://drop.io/mpfitdrop
>
>
>> Nice!
>
>
>> However, I see some points where the code could be improved for better
>
> reusability and maintainability.
>
>
>> If we lived in an ideal world with infinite time to polish everything,
>
> I'd like to see all of the points below addressed before landing this to
>
> Scipy. But since this would be lots of error-prone work, it's probably
>
> best to try to reach some compromise.
>
>
>> Given these constraints, I'd still like to see at least the coding style
>
> and error handling fixed (which probably are not too difficult to
>
> change), in addition to having better test coverage. The rest could come
>
> later, even if we accrue yet more technical debt with this...
>
>
>>
>> First, broad maintenance concerns:
>
>
>> - We already have `leastsq` from MINPACK. Having two MINPACK-derived
>
>  least squares fitting routines is not good.
>
>
>>  So, I'd perhaps like to see the `leastsq` core part extracted out of
>
>  MPFIT, and the MPFIT interface implemented on top of it as a thin
>
>  wrapper, or the other way around.
>
>
>>  Maybe, if the modifications made on MINPACK are small, they can be
>
>  backported to the Fortran code and MPFIT can be reimplemented on top
>
>  of `leastsq`.
>
>
>>  Any thoughts on this?
>
>
>> - What is the performance of the Python implementation as compared to the
>
>  Fortran code? For small data sets, the Python code is probably much
>
>  slower, but for large datasets is the performance is comparable?
>
>
>> - Much of the code is Fortran written in Python: long routines,
>
>  goto-equivalents, 6-letter variable names.
>
>
>>  Good commenting offset this, though.
>
>
>>
>> Then more specific points of refactoring:
>
>
>> - The code implements QR factorization with column pivoting from scratch.
>
>
>>  Perhaps some of this could be replaced with suitable LAPACK routines,
>
>  or with stuff from scipy.linalg. (Cf. DGEQPF, DGEQP3)
>
>
>>  I'm not sure whether there's something suitable for qrsolve in LAPACK;
>
>  the triangular system solver could be replaced with DTRTRS.
>
>
>>  Anyway, it might be useful to refactor qrfac and qrsolve out of MPFIT;
>
>  there may be other applications when it's useful to be able to solve
>
>  ||(A + I lambda) x - b||_2 = min! efficiently for multiple different
>
>  `lambda` in sequence.
>
>
>> - fdjac2 should probably be refactored out of MPFIT; other optimization
>
>  algorithms that need Jacobians computed by forward differences can then
>
>  use it. Do we already have something similar in Scipy already?
>
>
>> - `enorm` should be replaced with BLAS xNRM2; it does appropriate scaling
>
>  automatically and is probably much faster.
>
>
>>  enorm = scipy.lib.blas.get_blas_funcs(['nrm2'], [some_array])
>
>
>> - The long mpfit.__init__ routine should be split into smaller parts,
>
>  the function body is too long. I'm not sure exactly what parts, though.
>
>
>>  Perhaps at least the covariance matrix computation and input argument
>
>  parsing should be separated.
>
>
>> - "self.errmsg = 'ERROR ...; return'".
>
>  Probably should raise exceptions instead, at least for errors in input
>
>  parameters.
>
>
>>  In general, I think the error handling should make better use of the
>
>  Python exception and warning system; using `print` is not a correct
>
>  way to signal an error in library code.
>
>
>> - I'm not sure about implementing everything in a class. This tends to
>
>  couple tightly parts of the code that wouldn't really need such a strong
>
>  coupling.
>
>
>> - Does it work with complex-valued inputs?
>
>
>>
>> Numpy issues:
>
>
>> - Use numpy.finfo instead of machar
>
>
>> - Lots of nonzero/put/take combos, probably dates from Numeric.
>
>
>>  I think these should be replaced with boolean array indexing, to
>
>  enhance readability and performance.
>
>
>> - numpy.min/max -> amin/amax
>
>
>>
>> Minor stylistic issues in the code:
>
>
>> - `catch_msg` is not used for anything
>
>
>> - The original fortran documentation could be moved into the method
>
>  docstrings.
>
>
>> - "if foo:", not "if (foo):"
>
>
>> - if foo:
>
>     bar
>
>
>>  not
>
>
>>  if foo: bar
>
>
>> - "return foo", not "return(foo)"
>
>
>> - "# comment", not "## comment"
>
>
>> - It's not necessary to "import types"; you can just use
>
>  isinstance(x, float)
>
>
>>
>> [clip]
>
> > I have written a nose-test compatible test.  Could someone look at it
>
> > and tell me if it meets the scipy style before I continue adding more
>
> > tests from Craig's test-suite for the C version of his program?
>
>
>> The test style is OK; if it works with nosetests, it should be OK.
>
> Suggest using "import numpy as np", though, to stay in line with the rest
>
> of the code.
>
>
>> --
>
> Pauli Virtanen
>
>
>> _______________________________________________
>
> Scipy-dev mailing list
>
> Scipy-dev at scipy.org
>
> http://mail.scipy.org/mailman/listinfo/scipy-dev
>
>
> """
>
> We cleaned up most of these issues and would like SciPy to reconsider the
> package. Here is a short list of things we feel we accomplished:
>
> - Minor stylic issues cleaned up
> - Numpy issues (still a lot of nonzero calls but no more put/take combos)
> - There is an approx_jacobian that uses a forward diff method in
> scipy.optimize.fmin_slsq. We were not sure whether we should still factor
> out our jacobian approximation. We also just checked for jacobian
> approximators in SciPy - we haven't checked to see which is better, cleaner,
> faster, etc.
> - enorm fixed
> - long init routine split in half
> - Exceptions included (although a lot of them to show specific errors -
> this can be changed)
>
> Things that still need to be fixed:
>
> - Complex numbers. We didn't have this in mind in during implementation and
> still haven't thoroughly tested.
> - Still written as a class. What is the best alternative? Just a simple
> function?
> - QR from scratch.
> - Performance testing vs. Fortran code.
>
> Please let us know if there is anymore work you would like done on this
> optimizer. We use it for our own applications and it would be nice to have
> it included in SciPy for others to use. Thanks for your time.
>
> Bill Flynn
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20090819/af942acf/attachment.html>