Port features from SciPy's minpack version

[certik]

The version is here: https://github.com/scipy/scipy/tree/a0535bee9d4da58e6d41e963e907c1b9f7684585/scipy/optimize/minpack and we need to check to see what changes they have made to the original code and port those to this version.

Towards #14.

[ivan-pi]

I checked the history of that folder. The only major change is they added the recursive attribute to all subroutines to make sure they are re-entrant. The relevant issue is scipy/scipy#7999

If I understand correctly, by adding recursive the compiler must guarantee the subroutines are re-entrant. Quoting one the replies in that thread:

Many Fortran 77 compilears saved all local variables by default. gfortran still uses this misfeature, even for Fortran 90, albeit to a lesser extent and it can be controlled with compiler flags. Since the SciPy source cannot assume anything about the compiler and compiler flags passed to the Fortran compiler, we just have to assume that Fortran subroutines are never reentrant.

[milancurcic]

With F2018, procedures are recursive by default, nevertheless we should use recursive to ensure this on older compiler versions.

[ivan-pi]

Here's a list of things we could try to improve on based upon a search of the open issues containing the keyword "minpack" in the SciPy repository:

• Robust handling of NaN and Inf; scipy.optimize.leastsq: handling of inf and nan. scipy/scipy#4777, Ignore NaN with scipy.optimize.curve_fit scipy/scipy#11841
• Loss functions for robust least squares (e.g. filtering outliers); ENH: least_squares need 'huber' or other loss function in the method of 'lm' scipy/scipy#15162
• Count number of function evaluations: fsolve & root incorrect function-call count scipy/scipy#5369
• Better scaling of the objective function; related to fsolve can't find proper solution scipy/scipy#4227 and fmin_slsqp not converging to sensible solution scipy/scipy#7519
• Sparse Jacobians; Better routines for large-scale root finding (Trac #1790) scipy/scipy#2309, also see my suggestion in Specialization for banded Jacobians #1 about banded Jacobians. This is a bigger issue, because it means we also need specialized factorization routines for such matrices.
• Box constraints (I believe these are the same as lower and upper bounds on the solution); wishlist: root finding with box constraints scipy/scipy#5249

[certik]

@ivan-pi awesome! Let's implement those, and bring it to SciPy to use our newer improved version, which would fix a lot of their (reported) bugs.

[ivan-pi]

Perhaps more relevant than issues directly related to MINPACK is to look directly at the options available with other solvers in scipy.optimize.leastsq and scipy.optimize.root.

For least squares, the features missing in MINPACK compared to the other two SciPy solvers "trf" and "dogbox" are:

• Bound constraints
• Sparse Jacobian matrices
• Robust loss functions

The C++ library Ceres (BSD-licensed) implements all of these.

[xecej4]

Many Fortran 77 compilears saved all local variables by default. gfortran still uses this misfeature, even for Fortran 90, albeit to a lesser extent and it can be controlled with compiler flags.

That was stated in 2017. I do not recall Gfortran having that "misfeature", and recent versions of Gfortran do not do so. In fact, see this quote from the current Gfortran documentation:

"The default, which is -fautomatic, uses the stack for local variables smaller than the value given by -fmax-stack-var-size. Use the option -frecursive to use no static memory."

The Gfortran option "-fmax-stack-var-size=n" allows control over which arrays are put into the stack, and to avoid stack overflow. The parameter n has a default value of n = 65536 bytes. Thus, most local arrays of small or medium size are allocated on the stack without the need for any explicit request by the user.