fortran-lang / minpack

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

Provide Python bindings for Minpack #49

Closed awvwgk closed 2 years ago

awvwgk commented 2 years ago

Exports Python bindings via CFFI, allows installation of Python bindings with meson and setuptools.

Implementation

https://github.com/

Provided functions in minpack Python module

Follow-up PRs

Known issues

See #43

certik commented 2 years ago

Great job, thanks for pushing this!

I believe SciPy exclusively uses f2py for the wrappers? So I wonder if we should just use that, in order to be consistent.

@haozeke, what would you recommend?

But to get started, I think it doesn't matter how we do it, as long as it works, and then based on the feedback from the SciPy developers, we can adjust.

awvwgk commented 2 years ago

I believe SciPy exclusively uses f2py for the wrappers? So I wonder if we should just use that, in order to be consistent.

With f2py we would wrap the Fortran routines directly, rather than using the C bindings. The main drawback I see is that we can't pass-through data and exceptions via the Fortran routines easily (see #30).

certik commented 2 years ago

Yes, I agree we definitely want a robust C interface (header files) and we might as well use it for the Python wrappers.

I know @haozeke and I brainstormed this a few times, also in relation to have LFortran's pywrap generate the interface and if there is a way to join efforts with f2py. We can meet and brainstorm more.

In the meantime, I think doing what you are doing is fine. Let's get the wrappers working, and then we'll have something concrete in our hands, and experience with it, and we can see how to design a solution for long term.

Nicholaswogan commented 2 years ago

Are these binding compatible with numba compiled python? Can you do this:

import minpack.library
import numpy as np
from math import sqrt
import numba as nb

@nb.njit
def fcn(x, fvec) -> None:
    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)

@nb.njit
def test():
    res = minpack.library.hybrd1(fcn, x, fvec, tol)

test()

Compatibility with numba, or some python compiler, would be an improvement over the Scipy bindings.

awvwgk commented 2 years ago

Are these binding compatible with numba compiled python? Can you do this:

Seems like this is not possible right now, do you have an idea what I have to change for this to work?

❯ python test_numba.py
Traceback (most recent call last):
  File "/home/awvwgk/projects/src/git/minpack/test_numba.py", line 22, in <module>
    test()
  File "/home/awvwgk/software/opt/conda/envs/minpack/lib/python3.10/site-packages/numba/core/dispatcher.py", line 468, in _compile_for_args
    error_rewrite(e, 'typing')
  File "/home/awvwgk/software/opt/conda/envs/minpack/lib/python3.10/site-packages/numba/core/dispatcher.py", line 409, in error_rewrite
    raise e.with_traceback(None)
numba.core.errors.TypingError: Failed in nopython mode pipeline (step: nopython frontend)
Unknown attribute 'hybrd1' of type Module(<module 'minpack.library' from '/home/awvwgk/projects/src/git/minpack/python/minpack/library.py'>)

File "test_numba.py", line 20:
def test():
    res = minpack.library.hybrd1(fcn, x, fvec, tol)
    ^

During: typing of get attribute at /home/awvwgk/projects/src/git/minpack/test_numba.py (20)

File "test_numba.py", line 20:
def test():
    res = minpack.library.hybrd1(fcn, x, fvec, tol)
    ^
Nicholaswogan commented 2 years ago

Here is the relevant part of the numba docs: https://numba.pydata.org/numba-doc/latest/reference/pysupported.html?highlight=cffi

awvwgk commented 2 years ago

The Numba documentation is unfortunately not very useful due to its brevity. I did look a bit into the Numba source code trying to figure out what I'm supposed to do. I added a branch python-numba with some changes, but now I'm kind of stuck.

Current error status ``` ❯ python test_numba.py Traceback (most recent call last): File "/home/awvwgk/projects/src/git/minpack/test_numba.py", line 22, in test() File "/home/awvwgk/software/opt/conda/envs/minpack/lib/python3.10/site-packages/numba/core/dispatcher.py", line 468, in _compile_for_args error_rewrite(e, 'typing') File "/home/awvwgk/software/opt/conda/envs/minpack/lib/python3.10/site-packages/numba/core/dispatcher.py", line 409, in error_rewrite raise e.with_traceback(None) numba.core.errors.TypingError: Failed in nopython mode pipeline (step: nopython frontend) Failed in nopython mode pipeline (step: nopython frontend) Untyped global name 'minpack_hybrd1': Cannot determine Numba type of File "python/minpack/library.py", line 226: def hybrd1( minpack_hybrd1( ^ During: resolving callee type: type(CPUDispatcher()) During: typing of call at /home/awvwgk/projects/src/git/minpack/test_numba.py (20) During: resolving callee type: type(CPUDispatcher()) During: typing of call at /home/awvwgk/projects/src/git/minpack/test_numba.py (20) File "test_numba.py", line 20: def test(): res = minpack.library.hybrd1(fcn, x, fvec, tol) ^ ```

@Nicholaswogan if you have an idea what is needed to register this function with Numba let me know. Note that I decorated the function minpack_hybrd1 to get the necessary pre- and post-processing steps in place to pass-in the callback and pass-through the exceptions without disturbing the foreign runtimes.

awvwgk commented 2 years ago

Independent of the Numba compatibility this patch is now finished and ready for review. If we can figure the Numba compatibility out easily, I will add them here as well, but this doesn't affect the overall API design IMO.

Nicholaswogan commented 2 years ago

I'll play around with your branch and let you know if I figure out anything useful.

HaoZeke commented 2 years ago

I believe SciPy exclusively uses f2py for the wrappers? So I wonder if we should just use that, in order to be consistent.

With f2py we would wrap the Fortran routines directly, rather than using the C bindings. The main drawback I see is that we can't pass-through data and exceptions via the Fortran routines easily (see #30).

This should be possible in pyf files, but for now this approach seems like a good idea for a standalone library. Once these are in-place as @certik suggested, it would be good to check where the hand-generated bindings differ from the standard f2py generated signatures and determine from there which would be the best way forward for downstream projects like SciPy.

awvwgk commented 2 years ago

This should be possible in pyf files, but for now this approach seems like a good idea for a standalone library. Once these are in-place as @certik suggested, it would be good to check where the hand-generated bindings differ from the standard f2py generated signatures and determine from there which would be the best way forward for downstream projects like SciPy.

A problem no automatic wrapping can solve (at least the one I know of) is that the wrapped functions still look and feel like they are C or Fortran which is usually not really Pythonic.

I guess even f2py wrapping the routines directly would require an additional layer of Python to hide some of the implementation details which are mostly irrelevant in Python but important in a compiled language. I doubt that the resulting wrappers would be significantly smaller than what is in python/minpack/library.py already. I was actually surprised that one can produce effectively the same signature for the high and low level wrappers with some effort.

However, I prefer a route via C due to the general stability of the ABI in C, which is simply not available in Fortran. This becomes especially important when distributing (e.g. when packaging the Python bindings and the Fortran library separately).

The main issue in the current API I can see, might be that we have to modify the fvec / fjac in-place rather than return them. Not sure whether it is necessary to be conservative with memory usage here. For most of the tests I found that returning fvec and fjac requires longer callback implementations, but this is mainly due to the non-Pythonic for-loop usage and missing vectorization of the functions.

Nicholaswogan commented 2 years ago

I'm confused by how to test the python binding. I'm getting a segmentation fault.

I build the library:

meson builddir -Dpython=true
cd builddir  
meson compile

Then I copy pasted the built libminpack.dylib and _libminpack.cpython-39-darwin.so to the folder python/minpack/.

Then I added this directory to my runtime path

export DYLD_LIBRARY_PATH=/Users/nicholas/Downloads/minpack/python/minpack

The from the directory python/ I tried to run the following python script.

import pytest
import minpack.library
import numpy as np
from math import sqrt

def fcn(x, fvec) -> None:
    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)

assert minpack.library.hybrd1(fcn, x, fvec, tol) == 1

Segmentation fault occurs on the last line when minpack is called.

Nicholaswogan commented 2 years ago

Oh never mind. Its because i'm using a M1 mac, and minpack_capi.f90 uses nested functions, which don't work with current gfortran release on M1 mac. I won't be any help then.

certik commented 2 years ago

Segmentation fault occurs on the last line when minpack is called.

Unfortunately that is a common problem when the bindings are not compiled (ctypes, cffi?) unlike Cython or manual C/API. One has to go in, and debug it. One approach is to get a minimal example working on your machine, then slowly adding more complexity until you hit the bug. Probably some C signature incompatibility.

Nicholaswogan commented 2 years ago

@awvwgk, Current python API would take some work to be compatible with numba. Numba can not call any non-compiled python. So this sort of thing doesn't work, because test1 isn't compiled

@nb.njit
def test():
    test1()

def test1():
    print('Hi')

Looking at python/minpack/library.py, the callback scheme involves a lot of python.

I put together a wrapper hybrd and lmdif using ctypes a while back that is numba compatible: https://github.com/Nicholaswogan/NumbaMinpack . The API is a bit clunky, because I just couldn't figure out a better way.

awvwgk commented 2 years ago

Segmentation fault occurs on the last line when minpack is called.

Unfortunately that is a common problem when the bindings are not compiled (ctypes, cffi?) unlike Cython or manual C/API.

Even the C bindings should fail here, because GFortran for MacOS/Arm64 does not support nested function usage. If we decide on implementing #30 I can create an API which also works on MacOS/Arm64 before support in GFortran is implemented. Also, CFFI is compiled, just like Cython or building it manually in C.

Looking at python/minpack/library.py, the callback scheme involves a lot of python.

I put together a wrapper hybrd and lmdif using ctypes a while back that is numba compatible: https://github.com/Nicholaswogan/NumbaMinpack . The API is a bit clunky, because I just couldn't figure out a better way.

Indeed, I need some Python logic to ensure that we don't by accident throw an exception throw the Fortran runtime, as it would leave the stack in an undefined state. I have checked your ctypes based bindings already, we could use this as well if it is easier, however the logic on the Python side to handle exceptions would need to stay.

Edit: Looks like we are almost there, I now have to figure out how to properly register all CFFI procedures

Current error status ``` ❯ python test_numba.py Traceback (most recent call last): File "/home/awvwgk/projects/src/git/minpack/test_numba.py", line 1, in import minpack.library File "/home/awvwgk/projects/src/git/minpack/python/minpack/__init__.py", line 6, in from .library import ( File "/home/awvwgk/projects/src/git/minpack/python/minpack/library.py", line 203, in minpack_hybrd1 = cffi_callback(lib.func)(_libminpack.lib.minpack_hybrd1) File "/home/awvwgk/projects/src/git/minpack/python/minpack/library.py", line 174, in wrapper return dec(func, *args, **kwargs) File "/home/awvwgk/projects/src/git/minpack/python/minpack/library.py", line 199, in cffi_callback return njit(entry_point) File "/home/awvwgk/software/opt/conda/envs/minpack/lib/python3.10/site-packages/numba/core/decorators.py", line 258, in njit return jit(*args, **kws) File "/home/awvwgk/software/opt/conda/envs/minpack/lib/python3.10/site-packages/numba/core/decorators.py", line 179, in jit return wrapper(pyfunc) File "/home/awvwgk/software/opt/conda/envs/minpack/lib/python3.10/site-packages/numba/core/decorators.py", line 208, in wrapper disp = dispatcher(py_func=func, locals=locals, File "/home/awvwgk/software/opt/conda/envs/minpack/lib/python3.10/site-packages/numba/core/dispatcher.py", line 826, in __init__ pysig = utils.pysignature(py_func) File "/home/awvwgk/software/opt/conda/envs/minpack/lib/python3.10/inspect.py", line 3247, in signature return Signature.from_callable(obj, follow_wrapped=follow_wrapped, File "/home/awvwgk/software/opt/conda/envs/minpack/lib/python3.10/inspect.py", line 2995, in from_callable return _signature_from_callable(obj, sigcls=cls, File "/home/awvwgk/software/opt/conda/envs/minpack/lib/python3.10/inspect.py", line 2461, in _signature_from_callable return _signature_from_builtin(sigcls, obj, File "/home/awvwgk/software/opt/conda/envs/minpack/lib/python3.10/inspect.py", line 2271, in _signature_from_builtin raise ValueError("no signature found for builtin {!r}".format(func)) ValueError: no signature found for builtin ```
awvwgk commented 2 years ago

I updated the description with a bit more context on the implementation of the Python bindings. Let me know if there is anything you are not happy with, since we can refine the interface generation later, I hope nobody considers CFFI as a blocker.