search

fortran-lang / minpack

Modernized Minpack: for solving nonlinear equations and nonlinear least squares problems
https://fortran-lang.github.io/minpack/
Other
92 stars 19 forks source link

Goal: Get this minpack used by SciPy #14

Open certik opened 2 years ago

certik commented 2 years ago

SciPy contains minpack: https://github.com/scipy/scipy/tree/main/scipy/optimize/minpack; we should develop this fortran-lang/minpack in a way so that SciPy could use it. We should add some improvements (perhaps some new minimization algorithms), and then include them in SciPy. If we can pull this off, that would be a huge win, for both Fortran and SciPy. SciPy gets a maintainable version, we can all contribute to it in Fortran, and modern Fortran gets a "user" in a high profile Python library.

Fortran Forum discussion thread: https://fortran-lang.discourse.group/t/lets-get-fortran-lang-minpack-into-scipy/2722

jacobwilliams commented 2 years ago

One thing TODO: check to see what changes they have made to the original code and port those to this version. See #23.

milancurcic commented 2 years ago

I agree, this would be a huge boon for Fortran and SciPy. Great idea!

awvwgk commented 2 years ago

Here is a quick outline of the things I found in SciPy on this:

Note that there is a scipy.optimize.minpack2 module, this is interfaced via f2py

sakamoti commented 2 years ago

How about considering adding optional argument class(*) for fortran user?

In the f77 version, the user defined subroutine, which is the argument like hybrd1, is defined as func(n, x, fvec, iflag), etc. The user defined derived type cannot be passed. So, I use module variables or the internal subroutine to refer to the type.

By adding class(*) as an optional argument to func, for example, I think it is straightforward way to pass user-defined types. Personally, Such an additional argument class(*),optional will provide well userbility for fortran user.

Is it practically possible to add an optional argument like class(*),optional when considering the integration with Scipy?

jacobwilliams commented 2 years ago

See my comments in #4. I think that's the way to go, rather than having an optional argument.

awvwgk commented 2 years ago

There is also a relevant discussion on the callback design in https://fortran-lang.discourse.group/t/modern-fortran-interface-for-nlopt/2685/6 and https://github.com/grimme-lab/nlopt-f/issues/16. I think this is a quite fundamental issue, and I'm hitting against the limitation of the current callback when building the C API already.

sakamoti commented 2 years ago

@jacobwilliams Thank you for the link. I had missed that discussion. The code in the link is enough to serve the purpose.

@awvwgk I really appreciate it, because discussions like the one you linked to are hard to find outside of fortran-lang.

awvwgk commented 2 years ago

Is it practically possible to add an optional argument like class(*),optional when considering the integration with Scipy?

I would go so far to say that is necessary to pass additional data with the procedure pointer to allow integration with SciPy, Right now we have to rely on C and Fortran using similar enough calling conventions to use the same callback function interface. Once we switch to a more Fortran-friendly interface (assumed shape arrays instead of explicit shape arrays) this is no longer true.

Nicholaswogan commented 2 years ago

One other thing to add that I think is relevant. Currently released GNU compilers (e.g. GNU Fortran (Homebrew GCC 11.2.0_3) 11.2.0) do not support nested functions on Apple Silicon. This will change someday but it isn't clear when (see here). If Scipy relies on GNU compilers to build binaries, then passing data via nested functions will prevent building Apple Silicon binaries.

ilayn commented 2 years ago

Hi folks big fan of this initiative!

Let me know if you need any help from SciPy side please, @mdhaber would also be interested I presume, as the scipy.optimize resident and may involve others.

There is quite some legacy luggage with respect to xPACK routines in scipy. Hence if you think there is a more modern and "sane" way to wrap it, we are all ears. That is to say, please don't feel obliged to fit a circle to our strange square build mechanisms.

mdhaber commented 2 years ago

Yes, this sounds like it could be good. Have you found any issues that this would fix?

certik commented 2 years ago

@ilayn, @mdhaber thanks! Here is list of ideas that @ivan-pi collected that we could fix in minpack that would fix several SciPy issues:

If you have any other ideas, please let us know. We'll try to fix some of these, to "motivate" the transition to the new version, so that SciPy gains something concrete from it. Of those issues, which do you think would be the highest priority for SciPy?

mdhaber commented 2 years ago

Of those issues, which do you think would be the highest priority for SciPy?

Personally I'd suggest bug fixes first, then new features. Otherwise, I don't have a sense of priority, so I'd just go from easiest to hardest.

awvwgk commented 2 years ago

Is anyone opposed to adding meson build files to minpack? This way SciPy could easily give our minpack a try via a meson subproject now that the C API is merged.

certik commented 2 years ago

Yes, go ahead and add Meson!

ivan-pi commented 2 years ago

Given that SciPy moved to Meson it makes a lot of sense to support a build-system which makes it easier for them. 👍🏻

awvwgk commented 2 years ago

Sounds good, we now have meson support.

ilayn commented 2 years ago

@rgommers Do we need anything Meson specific with the modules? I think this is going to fit like a glove otherwise.

rgommers commented 2 years ago

No, things should work out of the box easily build-system wise. The work will be to do the integration step (differences in Fortran sources, f2py wrapper changes if needed, etc.).

Here is a quick outline of the things I found in SciPy on this: ....

Note that the header files, minpack.py, minpack2.py etc. are all private. The only thing we care about is the public functionality in the scipy.optimize namespace and that all the tests pass. Anything else is an implementation detail that can be changed as needed.

certik commented 2 years ago

Awesome, thanks @rgommers for the feedback. I think the move to Meson was a good move for SciPy, it should make this a lot easier.

awvwgk commented 2 years ago

I started working on Python bindings for Minpack in #49 using CFFI. To keep this thread on topic I opened a separate issue to discuss the general design of the Python API in #43, feel free to add links and resources to the description.

So far I only provide wrappers for the high-level functions (hybrd1, hybrj1, ...), with #49 the following snippet should work:

import numpy as np
from math import sqrt
from minpack import hybrd1

def fcn(x, fvec) -> None:
    """
    The problem is to determine the values of x(1), x(2), ..., x(9),
    which solve the system of tridiagonal equations.::

        (3-2*x(1))*x(1)           -2*x(2)                   = -1
                -x(i-1) + (3-2*x(i))*x(i)         -2*x(i+1) = -1, i=2-8
                                    -x(8) + (3-2*x(9))*x(9) = -1
    """

    for k in range(x.size):
        tmp = (3.0 - 2.0 * x[k]) * x[k]
        tmp1 = x[k - 1] if k > 0 else 0.0
        tmp2 = x[k + 1] if k < len(x) - 1 else 0.0
        fvec[k] = tmp - tmp1 - 2.0 * tmp2 + 1.0

x = np.array(9 * [-1.0])
fvec = np.zeros(9, dtype=np.float64)
tol = sqrt(np.finfo(np.float64).eps)

hybrd1(fcn, x, fvec, tol)

print(x)
# [-0.57065451 -0.68162834 -0.70173245
#  -0.70421294 -0.70136905 -0.69186564
#  -0.66579201 -0.5960342  -0.41641206]

Feedback is welcome. I have some experience in interop, but since I'm mostly working on the compiled language side, my Python skill is somewhat lacking behind ;).

awvwgk commented 2 years ago

Any thoughts on releasing the current version of minpack and start packaging it for conda-forge and other?

We could create a release candidate, but I doubt that there will be significant changes in the API in the next time.

certik commented 2 years ago

Sounds good to me.

xecej4 commented 2 years ago

After reading some of the bug reports on the use of Minpack in SciPy ( https://github.com/scipy/scipy/issues?page=2&q=minpack+bug ) , I came away with the feeling that bugs may be located in many places in the calling chain and the wrappers, or caused by compiling the Fortran part with incorrect compiler options. I suspect that most of the errors are not related to the stable Minpack Fortran code from Netlib, etc. However, if we replace the current Minpack sources with the new Fortran-Lang sources, and problems occur, the new Minpack code is likely to be suspected, even if the true culprit is the wrapper code or something else. Are there tests to verify the Fortran part independently, without the involvement of the wrapper code, Python, etc. ?

certik commented 2 years ago

@xecej4, good question. Yes, there are Fortran tests here: https://github.com/fortran-lang/minpack/tree/main/test and we can/will add more as needed. We also maintain Python wrappers in this repository here: https://github.com/fortran-lang/minpack/tree/main/python, we should add test for them as well. I suggest that we are responsible for the Fortran code and the wrappers and ensure everything works.

Then when a bug in SciPy happens, we simply reproduce it, test it, fix it. If things work in this repository, then we know the bug is caused by the way SciPy includes or compiles Minpack, either Fortran or the Python wrappers, so we can focus on just that and help them fix it.

That way there is a clean separation and it will workout great.

ilayn commented 2 years ago

Just to add some SciPy take, the first place we look at for faults are indeed our wrappers and other plumbing code that bridges over to compiled bits. Fortunately, there is not much in between and most bridges are quite battle tested. Thus if there is some trivial wrapping mistake we happen to find them quickly.

I, and maybe others would too, would like to help you folks with the tests and maybe carry them over to here from SciPy side to kickstart the testing suite.

Our biggest issue is that we keep track of patches of Fortran bits within SciPy with limited maintainer knowledge and not being able to judge fairly whether certain things are suitable for fixing-- say things like typed things with REAL*8 like stuff, threadsafe directives or other Fortran specific behaviors that we might not understand and/or evaluate accurately.

In that sense, while I don't foresee a big maintenance burden since it is stabilized by now for a long time already, there might be some bits that maybe modernized or replaced. It would be nice if we can follow your lead on that.

awvwgk commented 2 years ago

I just created a first release candidate at https://github.com/fortran-lang/minpack/releases/tag/v2.0.0-rc.1, including some Python wheels and sdists. Please let me know whether those work as expected.

arjenmarkus commented 2 years ago

Nice work, just a few remarks here: when I built the library and the test programs I got a number of warnings about array accesses out of bound from gfortran, While the message for "if ( k /= 1) temp1 = x(k - 1)" is clearly a false positive, the code for the test programs also contains constructions like:

do icase = 1, ncases+1

    if (icase == ncases+1) then
        write (nwrite, '(a,i3,a/)') '1SUMMARY OF ', lnp, ' TESTS OF

CHKDER' write (nwrite, '(a/)') ' NPROB N STATUS ERRMIN ERRMAX' do i = 1, lnp write (nwrite, '(i4, i6, 6x, l1, 3x, 2d15.7)') np(i), na(i), a(i), errmin(i), errmax(i) end do stop else nprob = nprobs(icase) n = ns(icase) ldfjac = n

Why do it that way? Why not simply put writing the summary AFTER the do-loop over the cases?

There is also a warning about the function dfloat(i): Warning: 'dfloat' declared at (1) may shadow the intrinsic of the same name. In order to call the intrinsic, explicit INTRINSIC declarations may be required. [-Wintrinsic-shadow] I am not sure if any of these would be of influence of influence on the acceptance of the code, but I thought I'd bring it up.

Op vr 20 mei 2022 om 21:53 schreef Sebastian Ehlert < @.***>:

I just created a first release candidate at https://github.com/fortran-lang/minpack/releases/tag/v2.0.0-rc.1, including some Python wheels and sdists. Please let me know whether those work as expected.

— Reply to this email directly, view it on GitHub https://github.com/fortran-lang/minpack/issues/14#issuecomment-1133277807, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAN6YR27O27IMO7FDIN4BYDVK7USXANCNFSM5NMJVQKQ . You are receiving this because you are subscribed to this thread.Message ID: @.***>

awvwgk commented 2 years ago

One thing where we could use some helping hands is the wheel builder, see .github/workflows/wheel.yml, which currently only creates Manylinux wheels (Linux/x86_64, CPython 3.7–3.10). Having MacOS and Windows wheels would be important for a PyPI upload once we release 2.0.0 (see #44). Building PyPy wheels on those platforms might be feasible as well.

rgommers commented 2 years ago

Maybe I'm missing some context @awvwgk, but I hope you are not working on wheels for SciPy's sake? We'd want to depend on this minpack as a git submodule, not as a separate dependency.

awvwgk commented 2 years ago

We'd want to depend on this minpack as a git submodule, not as a separate dependency.

A meson subproject seems preferable than. However, for channels like conda-forge I would prefer to not have SciPy vendor dependencies which can be packaged separately. Question is whether you want to use the provided Python bindings of this project or write your own bindings to the Minpack C-API.

rgommers commented 2 years ago

Sure, we can do a Meson subproject I think. That has been discussed for a few other git submodules as well - Boost in particular. Linux distros with strict rules may be happy with that. And I think we'd be happy to adopt the Python bindings from this project.

That said, I don't think we'd want a new runtime dependency for the scipy conda-forge package. But that's a thing that can be revisited later/separately.

h-vetinari commented 2 years ago

Hi, one of the conda-forge feedstock maintainers for scipy (not claiming to speak for the others).

A meson subproject seems preferable than. However, for channels like conda-forge I would prefer to not have SciPy vendor dependencies which can be packaged separately.

That's indeed the best practice for conda-forge.

Question is whether you want to use the provided Python bindings of this project or write your own bindings to the Minpack C-API.

I'm not 100% wure I understand the situation correctly, but from conda-forge side, it's not a problem to change build invocations and (less ideal but still possible) even patch things in or out, but I'd really like to keep the code-as-packaged as closely possible the same as what scipy tests and releases. IOW, switching the bindings should only happen if we can reasonably guarantee that they are using the exact same definition/release/commit