diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 4c4b63b..f0401d0 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -24,6 +24,14 @@ jobs: with: submodules: recursive + - name: Cache GFortran install + if: ${{ contains(matrix.os, 'windows') }} + id: cache + uses: actions/cache@v2 + with: + path: ./mingw-w64 + key: gcc-${{ matrix.gcc }}-${{ matrix.os }} + - name: Install GFortran (MacOS) if: ${{ contains(matrix.os, 'macos') }} run: | @@ -40,7 +48,7 @@ jobs: --slave /usr/bin/gcov gcov /usr/bin/gcov-${{ matrix.gcc }} - name: Install GFortran (Windows) - if: ${{ contains(matrix.os, 'windows') }} + if: ${{ contains(matrix.os, 'windows') && steps.cache.outputs.cache-hit != 'true' }} run: | Invoke-WebRequest -Uri ${{ env.DOWNLOAD }} -OutFile mingw-w64.zip Expand-Archive mingw-w64.zip @@ -70,6 +78,7 @@ jobs: if: ${{ matrix.build == 'meson' }} run: >- meson setup _build + --libdir=lib --prefix=${{ contains(matrix.os, 'windows') && '$pwd\_dist' || '$PWD/_dist' }} - name: Compile project (meson) @@ -84,6 +93,146 @@ jobs: if: ${{ matrix.build == 'meson' }} run: meson install -C _build --no-rebuild + - name: Create package (Unix) + if: ${{ matrix.build == 'meson' && ! contains(matrix.os, 'windows') }} + run: | + tar cvf ${{ env.OUTPUT }} _dist + xz -T0 ${{ env.OUTPUT }} + echo "MINPACK_OUTPUT=${{ env.OUTPUT }}.xz" >> $GITHUB_ENV + env: + OUTPUT: minpack-${{ matrix.os }}.tar + + - name: Create package (Windows) + if: ${{ matrix.build == 'meson' && contains(matrix.os, 'windows') }} + run: | + tar cvf ${{ env.OUTPUT }} _dist + xz -T0 ${{ env.OUTPUT }} + echo "MINPACK_OUTPUT=${{ env.OUTPUT }}.xz" | Out-File -FilePath $env:GITHUB_ENV -Encoding utf8 -Append + env: + OUTPUT: minpack-${{ matrix.os }}.tar + + - name: Upload package + if: ${{ matrix.build == 'meson' }} + uses: actions/upload-artifact@v2 + with: + name: ${{ env.MINPACK_OUTPUT }} + path: ${{ env.MINPACK_OUTPUT }} + + + Python: + needs: + - Build + runs-on: ${{ matrix.os }} + defaults: + run: + shell: ${{ contains(matrix.os, 'windows') && 'powershell' || 'bash -l {0}' }} + strategy: + fail-fast: false + matrix: + build: [meson] + os: [ubuntu-latest, macos-latest] + gcc: [10] + python: ['3.7', '3.8', '3.9'] + + # Additional test for setuptools build + include: + - build: setuptools + os: ubuntu-latest + gcc: 10 + python: '3.9' + + env: + FC: gfortran + CC: gcc + MINPACK_OUTPUT: minpack-${{ matrix.os }}.tar.xz + + steps: + - name: Checkout code + uses: actions/checkout@v2 + + - name: Cache GFortran install + if: ${{ contains(matrix.os, 'windows') }} + id: cache + uses: actions/cache@v2 + with: + path: ./mingw-w64 + key: gcc-${{ matrix.gcc }}-${{ matrix.os }} + + - name: Install dependencies + uses: mamba-org/provision-with-micromamba@main + with: + environment-file: config/ci/python-env.yaml + extra-specs: | + python=${{ matrix.python }} + + - name: Install GFortran (MacOS) + if: ${{ contains(matrix.os, 'macos') }} + run: | + brew install gcc@${{ matrix.gcc }} + ln -s /usr/local/bin/gfortran-${{ matrix.gcc }} /usr/local/bin/gfortran + ln -s /usr/local/bin/gcc-${{ matrix.gcc }} /usr/local/bin/gcc + + - name: Install GFortran (Linux) + if: ${{ contains(matrix.os, 'ubuntu') }} + run: | + sudo update-alternatives \ + --install /usr/bin/gcc gcc /usr/bin/gcc-${{ matrix.gcc }} 100 \ + --slave /usr/bin/gfortran gfortran /usr/bin/gfortran-${{ matrix.gcc }} \ + --slave /usr/bin/gcov gcov /usr/bin/gcov-${{ matrix.gcc }} + + - name: Install GFortran (Windows) + if: ${{ contains(matrix.os, 'windows') && steps.cache.outputs.cache-hit != 'true' }} + run: | + Invoke-WebRequest -Uri ${{ env.DOWNLOAD }} -OutFile mingw-w64.zip + Expand-Archive mingw-w64.zip + echo "$pwd\mingw-w64\mingw64\bin" | Out-File -FilePath $env:GITHUB_PATH -Encoding utf8 -Append + shell: pwsh + env: + DOWNLOAD: "https://github.com/brechtsanders/winlibs_mingw/releases/download/10.3.0-12.0.0-9.0.0-r2/winlibs-x86_64-posix-seh-gcc-10.3.0-mingw-w64-9.0.0-r2.zip" + + - name: Download package + uses: actions/download-artifact@v2 + with: + name: ${{ env.MINPACK_OUTPUT }} + + - name: Unpack package (Unix) + if: ${{ ! contains(matrix.os, 'windows') }} + run: | + tar xvf ${{ env.MINPACK_OUTPUT }} + echo "MINPACK_PREFIX=$PWD/_dist" >> $GITHUB_ENV + + - name: Unpack package (Windows) + if: ${{ contains(matrix.os, 'windows') }} + run: | + tar xvf ${{ env.MINPACK_OUTPUT }} + echo "MINPACK_OUTPUT=${{ env.OUTPUT }}.xz" | Out-File -FilePath $env:GITHUB_ENV -Encoding utf8 -Append + + - name: Install Python extension module (pip) + if: ${{ matrix.build == 'setuptools' }} + run: pip3 install . -vv + working-directory: python + env: + PKG_CONFIG_PATH: ${{ env.PKG_CONFIG_PATH }}:${{ env.MINPACK_PREFIX }}/lib/pkgconfig + LD_RUNPATH_SEARCH_PATH: ${{ env.MINPACK_PREFIX }}/lib + + - name: Install Python extension module (meson) + if: ${{ matrix.build == 'meson' }} + run: | + set -ex + meson setup _build --prefix=$CONDA_PREFIX --libdir=lib + meson compile -C _build + meson install -C _build + working-directory: python + env: + PKG_CONFIG_PATH: ${{ env.PKG_CONFIG_PATH }}:${{ env.MINPACK_PREFIX }}/lib/pkgconfig + + - name: Test Python API + run: pytest --pyargs minpack --cov=minpack -vv + env: + LD_LIBRARY_PATH: ${{ env.LD_LIBRARY_PATH }}:${{ env.MINPACK_PREFIX }}/lib + DYLD_LIBRARY_PATH: ${{ env.DYLD_LIBRARY_PATH }}:${{ env.MINPACK_PREFIX }}/lib + + Docs: runs-on: ubuntu-latest defaults: diff --git a/config/ci/python-env.yaml b/config/ci/python-env.yaml new file mode 100644 index 0000000..be9e59b --- /dev/null +++ b/config/ci/python-env.yaml @@ -0,0 +1,14 @@ +name: python +channels: + - conda-forge +dependencies: + - python + - pip + - pkgconfig + - pytest + - pytest-cov + - coverage + - cffi + - numpy + - meson + - ninja diff --git a/include/minpack.h b/include/minpack.h index f78a16b..91f817c 100644 --- a/include/minpack.h +++ b/include/minpack.h @@ -19,6 +19,16 @@ typedef void (*minpack_func)( int* /* iflag */, void* /* udata */); +#ifdef MINPACK_CFFI +extern "Python" void MINPACK_CALL +func( + int /* n */, + const double* /* x */, + double* /* fvec */, + int* /* iflag */, + void* /* udata */); +#endif + /* * the purpose of hybrd is to find a zero of a system of * n nonlinear functions in n variables by a modification @@ -84,6 +94,18 @@ typedef void (*minpack_fcn_hybrj)( int* /* iflag */, void* /* udata */); +#ifdef MINPACK_CFFI +extern "Python" void MINPACK_CALL +fcn_hybrj( + int /* n */, + const double* /* x */, + double* /* fvec */, + double* /* fjac */, + int /* ldfjac */, + int* /* iflag */, + void* /* udata */); +#endif + /* * the purpose of hybrj is to find a zero of a system of * n nonlinear functions in n variables by a modification @@ -148,6 +170,19 @@ typedef void (*minpack_fcn_lmder)( int* /* iflag */, void* /* udata */); +#ifdef MINPACK_CFFI +extern "Python" void MINPACK_CALL +fcn_lmder( + int /* m */, + int /* n */, + const double* /* x */, + double* /* fvec */, + double* /* fjac */, + int /* ldfjac */, + int* /* iflag */, + void* /* udata */); +#endif + /* * the purpose of lmder is to minimize the sum of the squares of * m nonlinear functions in n variables by a modification of @@ -213,6 +248,17 @@ typedef void (*minpack_func2)( int* /* iflag */, void* /* udata */); +#ifdef MINPACK_CFFI +extern "Python" void MINPACK_CALL +func2( + int /* m */, + int /* n */, + const double* /* x */, + double* /* fvec */, + int* /* iflag */, + void* /* udata */); +#endif + /* * the purpose of lmdif is to minimize the sum of the squares of * m nonlinear functions in n variables by a modification of @@ -279,6 +325,18 @@ typedef void (*minpack_fcn_lmstr)( int* /* iflag */, void* /* udata */); +#ifdef MINPACK_CFFI +extern "Python" void MINPACK_CALL +fcn_lmstr( + int /* m */, + int /* n */, + const double* /* x */, + double* /* fvec */, + double* /* fjrow */, + int* /* iflag */, + void* /* udata */); +#endif + /* * the purpose of lmstr is to minimize the sum of the squares of * m nonlinear functions in n variables by a modification of diff --git a/meson.build b/meson.build index 52a600d..94eb4d0 100644 --- a/meson.build +++ b/meson.build @@ -8,6 +8,7 @@ project( 'buildtype=debugoptimized', ], ) +has_cc = add_languages('c', required: get_option('python'), native: false) minpack_lib = library( meson.project_name(), @@ -59,3 +60,7 @@ subdir('examples') # add the testsuite subdir('test') + +if get_option('python') + subdir('python/minpack') +endif diff --git a/meson_options.txt b/meson_options.txt new file mode 100644 index 0000000..df0f2b1 --- /dev/null +++ b/meson_options.txt @@ -0,0 +1,12 @@ +option( + 'python', + type: 'boolean', + value: false, + description: 'Build Python extension module', +) +option( + 'python_version', + type: 'string', + value: 'python3', + description: 'Python version to link against.', +) diff --git a/python/.gitignore b/python/.gitignore new file mode 100644 index 0000000..9842ebd --- /dev/null +++ b/python/.gitignore @@ -0,0 +1,172 @@ +# Prerequisites +*.d + +# Compiled Object files +*.slo +*.lo +*.o +*.obj + +# Precompiled Headers +*.gch +*.pch + +# Compiled Dynamic libraries +*.so +*.dylib +*.dll + +# Fortran module files +*.mod +*.smod + +# Compiled Static libraries +*.lai +*.la +*.a +*.lib + +# Executables +*.exe +*.out +*.app + +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] +*$py.class + +# Distribution / packaging +.Python +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +wheels/ +share/python-wheels/ +*.egg-info/ +.installed.cfg +*.egg +MANIFEST + +# PyInstaller +# Usually these files are written by a python script from a template +# before PyInstaller builds the exe, so as to inject date/other infos into it. +*.manifest +*.spec + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt + +# Unit test / coverage reports +htmlcov/ +.tox/ +.nox/ +.coverage +.coverage.* +.cache +nosetests.xml +coverage.xml +*.cover +*.py,cover +.hypothesis/ +.pytest_cache/ +cover/ + +# Translations +*.mo +*.pot + +# Django stuff: +*.log +local_settings.py +db.sqlite3 +db.sqlite3-journal + +# Flask stuff: +instance/ +.webassets-cache + +# Scrapy stuff: +.scrapy + +# Sphinx documentation +docs/_build/ + +# PyBuilder +.pybuilder/ +target/ + +# Jupyter Notebook +.ipynb_checkpoints + +# IPython +profile_default/ +ipython_config.py + +# pyenv +# For a library or package, you might want to ignore these files since the code is +# intended to run in multiple environments; otherwise, check them in: +# .python-version + +# pipenv +# According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control. +# However, in case of collaboration, if having platform-specific dependencies or dependencies +# having no cross-platform support, pipenv may install dependencies that don't work, or not +# install all needed dependencies. +#Pipfile.lock + +# PEP 582; used by e.g. github.com/David-OConnor/pyflow +__pypackages__/ + +# Celery stuff +celerybeat-schedule +celerybeat.pid + +# SageMath parsed files +*.sage.py + +# Environments +.env +.venv +env/ +venv/ +ENV/ +env.bak/ +venv.bak/ + +# Spyder project settings +.spyderproject +.spyproject + +# Rope project settings +.ropeproject + +# mkdocs documentation +/site + +# mypy +.mypy_cache/ +.dmypy.json +dmypy.json + +# Pyre type checker +.pyre/ + +# pytype static type analyzer +.pytype/ + +# Cython debug symbols +cython_debug/ + +# Directories +/build*/ +/_*/ diff --git a/python/README.rst b/python/README.rst new file mode 100644 index 0000000..dc9e38c --- /dev/null +++ b/python/README.rst @@ -0,0 +1,85 @@ +Minpack Python bindings +======================= + +Python bindings for Minpack. + + +Building the extension module +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +To perform a build some version of ``minpack`` has to be available on your system and preferably findable by ``pkg-config``. +Try to find a ``minpack`` installation you build against first with + +.. code:: sh + + pkg-config --modversion minpack + +Adjust the ``PKG_CONFIG_PATH`` environment variable to include the correct directories to find the installation if necessary. + + +Using pip +^^^^^^^^^ + +This project support installation with pip as an easy way to build the Python API. + +- C compiler to build the C-API and compile the extension module (the compiler name should be exported in the ``CC`` environment variable) +- Python 3.6 or newer +- The following Python packages are required additionally + + - `cffi `_ + - `numpy `_ + - `pkgconfig `_ (setup only) + +Make sure to have your C compiler set to the ``CC`` environment variable + +.. code:: sh + + export CC=gcc + +Install the project with pip + +.. code:: sh + + pip install . + +If you already have a ``minpack`` installation, *e.g.* from conda-forge, you can build the Python extension module directly without cloning this repository + +.. code:: sh + + pip install "https://github.com/fortran-lang/minpack/archive/refs/heads/main.zip#egg=minpack&subdirectory=python" + + + +Using meson +^^^^^^^^^^^ + +This directory contains a separate meson build file to allow the out-of-tree build of the CFFI extension module. +The out-of-tree build requires + +- C compiler to build the C-API and compile the extension module +- `meson `_ version 0.53 or newer +- a build-system backend, *i.e.* `ninja `_ version 1.7 or newer +- Python 3.6 or newer with the `CFFI `_ package installed + +Setup a build with + +.. code:: sh + + meson setup _build -Dpython_version=$(which python3) + +The Python version can be used to select a different Python version, it defaults to ``'python3'``. +Python 2 is not supported with this project, the Python version key is meant to select between several local Python 3 versions. + +Compile the project with + +.. code:: sh + + meson compile -C _build + +The extension module is now available in ``_build/minpack/_libminpack.*.so``. +You can install as usual with + +.. code:: sh + + meson configure _build --prefix=/path/to/install + meson install -C _build diff --git a/python/ffi-builder.py b/python/ffi-builder.py new file mode 100644 index 0000000..d221f31 --- /dev/null +++ b/python/ffi-builder.py @@ -0,0 +1,83 @@ +""" +FFI builder module for minpack for usage from meson and from setup.py. + +Since meson has the full knowledge about the build, it will handle +the generation of the C definitions in the meson.build file rather +than in the FFI builder. This allows to correctly keep track of +dependencies and updates in the build process. + +For setup.py we have to do the preprocessing ourselves here, this +requires us to use the C compiler to preprocess the header file +of minpack because the CFFI C parser cannot handle certain C +preprocessor constructs. Also, we cannot rely on an external build +system fixing dependencies for us and therefore have to find those +ourselves using pkg-config. +""" + +import os +import cffi + +library = "minpack" +include_header = '#include "minpack.h"' +prefix_var = "MINPACK_PREFIX" +if prefix_var not in os.environ: + prefix_var = "CONDA_PREFIX" + +if __name__ == "__main__": + import sys + + kwargs = dict(libraries=[library]) + + header_file = sys.argv[1] + module_name = sys.argv[2] + + with open(header_file) as f: + cdefs = f.read() +else: + import subprocess + + try: + import pkgconfig + + if not pkgconfig.exists(library): + raise ModuleNotFoundError("Unable to find pkg-config package 'minpack'") + if pkgconfig.installed(library, "< 2.0.0"): + raise Exception( + "Installed 'minpack' version is too old, 2.0.0 or newer is required" + ) + + kwargs = pkgconfig.parse(library) + cflags = pkgconfig.cflags(library).split() + + except ModuleNotFoundError: + kwargs = dict(libraries=[library]) + cflags = [] + if prefix_var in os.environ: + prefix = os.environ[prefix_var] + kwargs.update( + include_dirs=[os.path.join(prefix, "include")], + library_dirs=[os.path.join(prefix, "lib")], + runtime_library_dirs=[os.path.join(prefix, "lib")], + ) + cflags.append("-I" + os.path.join(prefix, "include")) + + cc = os.environ["CC"] if "CC" in os.environ else "cc" + + module_name = "minpack._libminpack" + + p = subprocess.Popen( + [cc, *cflags, "-DMINPACK_CFFI=1", "-E", "-"], + stdin=subprocess.PIPE, + stdout=subprocess.PIPE, + stderr=subprocess.PIPE, + ) + out, err = p.communicate(include_header.encode()) + + cdefs = out.decode() + +ffibuilder = cffi.FFI() +ffibuilder.set_source(module_name, include_header, **kwargs) +ffibuilder.cdef(cdefs) + +if __name__ == "__main__": + ffibuilder.distutils_extension(".") diff --git a/python/meson.build b/python/meson.build new file mode 100644 index 0000000..91d3eaa --- /dev/null +++ b/python/meson.build @@ -0,0 +1,14 @@ +project( + 'minpack', + 'c', + meson_version: '>=0.53', + default_options: [ + 'buildtype=debugoptimized', + ], +) +install = true + +minpack_dep = dependency('minpack', version: '>=2.0.0') +minpack_header = files('../include/minpack.h') + +subdir('minpack') diff --git a/python/meson_options.txt b/python/meson_options.txt new file mode 100644 index 0000000..7c91d21 --- /dev/null +++ b/python/meson_options.txt @@ -0,0 +1,6 @@ +option( + 'python_version', + type: 'string', + value: 'python3', + description: 'Python version to link against.', +) diff --git a/python/minpack/__init__.py b/python/minpack/__init__.py new file mode 100644 index 0000000..bd87960 --- /dev/null +++ b/python/minpack/__init__.py @@ -0,0 +1,12 @@ +# -*- coding: utf-8 -*- +""" +Python bindings for the Minpack library. +""" + +from .library import ( + hybrd1, + hybrj1, + lmdif1, + lmder1, + lmstr1, +) diff --git a/python/minpack/library.py b/python/minpack/library.py new file mode 100644 index 0000000..702211d --- /dev/null +++ b/python/minpack/library.py @@ -0,0 +1,370 @@ +# -*- coding: utf-8 -*- +""" +Low-level Python bindings to the Minpack library. + +This module forwards the CFFI generated bindings to the Minpack library and provides +a Pythonic interface to the C API. +""" + +import numpy as np +import functools + +from ._libminpack import ffi, lib + + +class MinpackInputError(Exception): + """ + Exception raised when Minpack input is invalid. + """ + pass + + +class UserData: + """ + Carry Python callable object through callback and propagate exceptions without + disturbing foreign runtime. + """ + + def __init__(self, fcn): + self.fcn = fcn + self.exception = None + + +@ffi.def_extern() +def func(n, x, fvec, iflag, data) -> None: + """ + Entry point for callback from minpack_hybrd and minpack_hybrd1 library functions. + Restores type information for NDArray objects and calls the user-provided callback. + """ + if iflag[0] <= 0: return + handle: UserData = ffi.from_handle(data) + try: + handle.fcn( + np.frombuffer(ffi.buffer(x, n * real.itemsize), dtype=real), + np.frombuffer(ffi.buffer(fvec, n * real.itemsize), dtype=real), + ) + except BaseException as e: + iflag[0] = -1 + handle.exception = e + + +@ffi.def_extern() +def fcn_hybrj(n, x, fvec, fjac, ldfjac, iflag, data) -> None: + """ + Entry point for callback from minpack_hybrj and minpack_hybrj1 library functions. + Restores type information for NDArray objects and calls the user-provided callback. + """ + if iflag[0] <= 0: return + handle: UserData = ffi.from_handle(data) + try: + fjac = np.frombuffer(ffi.buffer(fjac, ldfjac * n * real.itemsize), dtype=real) + handle.fcn( + np.frombuffer(ffi.buffer(x, n * real.itemsize), dtype=real), + np.frombuffer(ffi.buffer(fvec, n * real.itemsize), dtype=real), + np.reshape(fjac, (n, ldfjac)), + iflag[0] == 2, + ) + except BaseException as e: + iflag[0] = -1 + handle.exception = e + + +@ffi.def_extern() +def fcn_lmder(m, n, x, fvec, fjac, ldfjac, iflag, data) -> None: + """ + Entry point for callback from minpack_lmder and minpack_lmder1 library functions. + Restores type information for NDArray objects and calls the user-provided callback. + """ + if iflag[0] <= 0: return + handle: UserData = ffi.from_handle(data) + try: + fjac = np.frombuffer(ffi.buffer(fjac, ldfjac * n * real.itemsize), dtype=real) + handle.fcn( + np.frombuffer(ffi.buffer(x, n * real.itemsize), dtype=real), + np.frombuffer(ffi.buffer(fvec, m * real.itemsize), dtype=real), + np.reshape(fjac, (n, ldfjac)), + iflag[0] == 2, + ) + except BaseException as e: + iflag[0] = -1 + handle.exception = e + + +@ffi.def_extern() +def func2(m, n, x, fvec, iflag, data) -> None: + """ + Entry point for callback from minpack_lmdif and minpack_lmdif1 library functions. + Restores type information for NDArray objects and calls the user-provided callback. + """ + if iflag[0] <= 0: return + handle: UserData = ffi.from_handle(data) + try: + handle.fcn( + np.frombuffer(ffi.buffer(x, n * real.itemsize), dtype=real), + np.frombuffer(ffi.buffer(fvec, m * real.itemsize), dtype=real), + ) + except BaseException as e: + iflag[0] = -1 + handle.exception = e + + +@ffi.def_extern() +def fcn_lmstr(m, n, x, fvec, fjrow, ldfjac, iflag, data) -> None: + """ + Entry point for callback from minpack_lmstr and minpack_lmstr1 library functions. + Restores type information for NDArray objects and calls the user-provided callback. + """ + if iflag[0] <= 0: return + handle: UserData = ffi.from_handle(data) + try: + handle.fcn( + np.frombuffer(ffi.buffer(x, n * real.itemsize), dtype=real), + np.frombuffer(ffi.buffer(fvec, m * real.itemsize), dtype=real), + np.frombuffer(ffi.buffer(fjrow, n * real.itemsize), dtype=real), + iflag[0] - 2 if iflag[0] > 1 else None, + ) + except BaseException as e: + iflag[0] = -1 + handle.exception = e + + +def extern_python(dec): + """ + Meta-decorator to attach a CFFI extern "Python" callback to a decorator + handling the Python side of the callback. + """ + + def layer(*args, **kwargs): + def wrapper(func): + return dec(func, *args, **kwargs) + + return wrapper + + return layer + + +@extern_python +def cffi_callback(func, callback): + """ + Attach Python callback to a library function with extern "Python" callback. + + This decorator wraps the user-provided Python callback in a `UserData` object + to carry it through the foreign runtime. It also propagates exceptions from + the user-provided Python callback back through the foreign runtime and re-raises. + """ + + @functools.wraps(func) + def entry_point(fcn, *args): + data = UserData(fcn) + handle = ffi.new_handle(data) + func(callback, *args, handle) + if data.exception is not None: + raise data.exception + + return entry_point + + +real = np.dtype("f8") +minpack_hybrd1 = cffi_callback(lib.func)(lib.minpack_hybrd1) +minpack_hybrd = cffi_callback(lib.func)(lib.minpack_hybrd) +minpack_hybrj1 = cffi_callback(lib.fcn_hybrj)(lib.minpack_hybrj1) +minpack_hybrj = cffi_callback(lib.fcn_hybrj)(lib.minpack_hybrj) +minpack_lmder1 = cffi_callback(lib.fcn_lmder)(lib.minpack_lmder1) +minpack_lmder = cffi_callback(lib.fcn_lmder)(lib.minpack_lmder) +minpack_lmdif1 = cffi_callback(lib.func2)(lib.minpack_lmdif1) +minpack_lmdif = cffi_callback(lib.func2)(lib.minpack_lmdif) +minpack_lmstr1 = cffi_callback(lib.fcn_lmstr)(lib.minpack_lmstr1) +minpack_lmstr = cffi_callback(lib.fcn_lmstr)(lib.minpack_lmstr) +minpack_chkder = lib.minpack_chkder + + +def hybrd1(fcn, x, fvec, tol) -> int: + """ + Find a zero of a system of n nonlinear functions in n variables + by a modification of the Powell hybrid method. + This is done by using the more general nonlinear equation solver `hybrd`. + The user must provide a subroutine which calculates the functions. + The Jacobian is then calculated by a forward-difference approximation. + """ + n = x.size + lwa = n * (3 * n + 13) // 2 + wa = np.zeros(lwa, dtype=real) + info = ffi.new("int *") + + minpack_hybrd1( + fcn, + n, + ffi.cast("double*", x.ctypes.data), + ffi.cast("double*", fvec.ctypes.data), + tol, + info, + ffi.cast("double*", wa.ctypes.data), + lwa, + ) + if info[0] == 0: + raise MinpackInputError("Improper input parameters in hybrd1") + return info[0] + + +def hybrj1(fcn, x, fvec, fjac, tol) -> int: + """ + Find a zero of a system of n nonlinear functions in n variables + by a modification of the Powell hybrid method. + This is done by using the more general nonlinear equation solver `hybrj`. + The user must provide a subroutine which calculates the functions + and the Jacobian. + """ + n = x.size + lwa = (n * (n + 13)) // 2 + wa = np.zeros(lwa, dtype=real) + info = ffi.new("int *") + + minpack_hybrj1( + fcn, + n, + ffi.cast("double*", x.ctypes.data), + ffi.cast("double*", fvec.ctypes.data), + ffi.cast("double*", fjac.ctypes.data), + n, + tol, + info, + ffi.cast("double*", wa.ctypes.data), + lwa, + ) + if info[0] == 0: + raise MinpackInputError("Improper input parameters in hybrj1") + return info[0] + + +def lmder1(fcn, x, fvec, fjac, tol) -> int: + """ + Minimize the sum of the squares of m nonlinear functions in n variables + by a modification of the Levenberg-Marquardt algorithm. + This is done by using the more general least-squares solver `lmder`. + The user must provide a subroutine which calculates the functions and the Jacobian. + """ + n = x.size + m = fvec.size + lwa = 5 * n + m + wa = np.zeros(lwa, dtype=real) + ipvt = np.zeros(n, dtype=np.int32) + info = ffi.new("int *") + + minpack_lmder1( + fcn, + m, + n, + ffi.cast("double*", x.ctypes.data), + ffi.cast("double*", fvec.ctypes.data), + ffi.cast("double*", fjac.ctypes.data), + m, + tol, + info, + ffi.cast("int*", ipvt.ctypes.data), + ffi.cast("double*", wa.ctypes.data), + lwa, + ) + if info[0] == 0: + raise MinpackInputError("Improper input parameters in lmder1") + return info[0] + + +def lmdif1(fcn, x, fvec, tol) -> int: + """ + Minimize the sum of the squares of m nonlinear functions in n variables + by a modification of the Levenberg-Marquardt algorithm. + This is done by using the more general least-squares solver `lmdif`. + The user must provide a subroutine which calculates the functions. + The jacobian is then calculated by a forward-difference approximation. + """ + n = x.size + m = fvec.size + lwa = m * n + 5 * n + m + wa = np.zeros(lwa, dtype=real) + ipvt = np.zeros(n, dtype=np.int32) + info = ffi.new("int *") + + minpack_lmdif1( + fcn, + m, + n, + ffi.cast("double*", x.ctypes.data), + ffi.cast("double*", fvec.ctypes.data), + tol, + info, + ffi.cast("int*", ipvt.ctypes.data), + ffi.cast("double*", wa.ctypes.data), + lwa, + ) + if info[0] == 0: + raise MinpackInputError("Improper input parameters in lmdif1") + return info[0] + + +def lmstr1(fcn, x, fvec, fjac, tol) -> int: + """ + Minimize the sum of the squares of m nonlinear functions in n variables by + a modification of the Levenberg-Marquardt algorithm which uses minimal storage. + This is done by using the more general least-squares solver `lmstr`. + The user must provide a subroutine which calculates the functions and + the rows of the Jacobian. + """ + n = x.size + m = fvec.size + lwa = 5 * n + m + wa = np.zeros(lwa, dtype=real) + ipvt = np.zeros(n, dtype=np.int32) + info = ffi.new("int *") + + minpack_lmstr1( + fcn, + m, + n, + ffi.cast("double*", x.ctypes.data), + ffi.cast("double*", fvec.ctypes.data), + ffi.cast("double*", fjac.ctypes.data), + n, + tol, + info, + ffi.cast("int*", ipvt.ctypes.data), + ffi.cast("double*", wa.ctypes.data), + lwa, + ) + if info[0] == 0: + raise MinpackInputError("Improper input parameters in lmstr1") + return info[0] + + +def chkder(x, fvec, fjac, xp, fvecp, check, error): + """ + This subroutine checks the gradients of m nonlinear functions + in n variables, evaluated at a point x, for consistency with + the functions themselves. + + The subroutine does not perform reliably if cancellation or + rounding errors cause a severe loss of significance in the + evaluation of a function. Therefore, none of the components + of x should be unusually small (in particular, zero) or any + other value which may cause loss of significance. + """ + if not fvec.size == fjac.shape[-1] == fvecp.size == error.size: + raise ValueError("fvec, fjac, fvecp, error must have the same size") + if not x.size == fjac.shape[0] == xp.size: + raise ValueError("x, fjac, xp must have the same size") + + m = fvec.size + n = x.size + ldfjac = fjac.shape[-1] + + minpack_chkder( + m, + n, + ffi.cast("double*", x.ctypes.data), + ffi.cast("double*", fvec.ctypes.data), + ffi.cast("double*", fjac.ctypes.data), + ldfjac, + ffi.cast("double*", xp.ctypes.data), + ffi.cast("double*", fvecp.ctypes.data), + 2 if check else 1, + ffi.cast("double*", error.ctypes.data), + ) diff --git a/python/minpack/meson.build b/python/minpack/meson.build new file mode 100644 index 0000000..f516e42 --- /dev/null +++ b/python/minpack/meson.build @@ -0,0 +1,43 @@ +cc = meson.get_compiler('c') +ext_module = '_libminpack' + +pymod = import('python') +python = pymod.find_installation( + get_option('python_version'), + modules: [ + 'cffi', + ], +) +python_dep = python.dependency(required: true) + +# Python's CFFI is horrible in working with preprocessor statements, +# therefore, we have to preprocess the header before passing it to the ffibuilder +minpack_pp = configure_file( + command: [cc, '-DMINPACK_CFFI=1', '-E', '@INPUT@'], + input: minpack_header[0], + output: '@0@.h'.format(ext_module), + capture: true, +) + +# This is the actual out-of-line API processing of the ffibuilder +minpack_cffi_srcs = configure_file( + command: [python, files('..'/'ffi-builder.py'), '@INPUT@', '@BASENAME@'], + input: minpack_pp, + output: '@BASENAME@.c', +) + +# Actual generation of the Python extension +minpack_pyext = python.extension_module( + ext_module, + minpack_cffi_srcs, + dependencies: [minpack_dep, python_dep], + install: true, + subdir: 'minpack', +) + +python.install_sources( + '__init__.py', + 'library.py', + 'test_library.py', + subdir: 'minpack', +) diff --git a/python/minpack/test_library.py b/python/minpack/test_library.py new file mode 100644 index 0000000..0a49fee --- /dev/null +++ b/python/minpack/test_library.py @@ -0,0 +1,229 @@ +# -*- coding: utf-8 -*- + +import pytest +import minpack.library +import numpy as np +from math import sqrt + + +def test_hybrd1(): + 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 + + assert pytest.approx(x, abs=10 * tol) == [ + -0.5706545, + -0.6816283, + -0.7017325, + -0.7042129, + -0.7013690, + -0.6918656, + -0.6657920, + -0.5960342, + -0.4164121, + ] + + +def test_hybrd1_exception(): + + class DummyException(Exception): ... + + def fcn(x, fvec) -> None: + raise DummyException() + + x = np.array(9 * [-1.0]) + fvec = np.zeros(9, dtype=np.float64) + tol = sqrt(np.finfo(np.float64).eps) + + with pytest.raises(DummyException): + minpack.library.hybrd1(fcn, x, fvec, tol) + + +def test_hybrj1(): + def fcn(x, fvec, fjac, jacobian: bool) -> None: + + if jacobian: + for k in range(x.size): + for j in range(x.size): + fjac[k, j] = 0.0 + fjac[k, k] = 3.0 - 4.0 * x[k] + if k > 0: + fjac[k, k - 1] = -1.0 + if k < x.size - 1: + fjac[k, k + 1] = -2.0 + else: + 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) + fjac = np.zeros((9, 9), dtype=np.float64) + tol = sqrt(np.finfo(np.float64).eps) + + assert minpack.library.hybrj1(fcn, x, fvec, fjac, tol) == 1 + + assert pytest.approx(x, abs=10 * tol) == [ + -0.5706545, + -0.6816283, + -0.7017325, + -0.7042129, + -0.7013690, + -0.6918656, + -0.6657920, + -0.5960342, + -0.4164121, + ] + + +def test_hybrj1_exception(): + + class DummyException(Exception): ... + + def fcn(x, fvec, fjac, jacobian) -> None: + raise DummyException() + + x = np.array(9 * [-1.0]) + fvec = np.zeros(9, dtype=np.float64) + fjac = np.zeros((9, 9), dtype=np.float64) + tol = sqrt(np.finfo(np.float64).eps) + + with pytest.raises(DummyException): + minpack.library.hybrj1(fcn, x, fvec, fjac, tol) + + +def test_lmder1(): + y = np.array( + [ + 1.4e-1, + 1.8e-1, + 2.2e-1, + 2.5e-1, + 2.9e-1, + 3.2e-1, + 3.5e-1, + 3.9e-1, + 3.7e-1, + 5.8e-1, + 7.3e-1, + 9.6e-1, + 1.34e0, + 2.1e0, + 4.39e0, + ] + ) + + def fcn(x, fvec, fjac, jacobian: bool) -> None: + if jacobian: + for i in range(fvec.size): + tmp1, tmp2 = i + 1, 16 - i - 1 + tmp3 = tmp2 if i >= 8 else tmp1 + tmp4 = (x[1] * tmp2 + x[2] * tmp3) ** 2 + fjac[0, i] = -1.0 + fjac[1, i] = tmp1 * tmp2 / tmp4 + fjac[2, i] = tmp1 * tmp3 / tmp4 + else: + for i in range(fvec.size): + tmp1, tmp2 = i + 1, 16 - i - 1 + tmp3 = tmp2 if i >= 8 else tmp1 + fvec[i] = y[i] - (x[0] + tmp1 / (x[1] * tmp2 + x[2] * tmp3)) + + x = np.array([1.0, 1.0, 1.0]) + fvec = np.zeros(15, dtype=np.float64) + fjac = np.zeros((3, 15), dtype=np.float64) + tol = sqrt(np.finfo(np.float64).eps) + + xp = np.zeros(3, dtype=np.float64) + fvecp = np.zeros(15, dtype=np.float64) + err = np.zeros(15, dtype=np.float64) + minpack.library.chkder(x, fvec, fjac, xp, fvecp, False, err) + fcn(x, fvec, fjac, False) + fcn(x, fvec, fjac, True) + fcn(xp, fvecp, fjac, False) + minpack.library.chkder(x, fvec, fjac, xp, fvecp, True, err) + + assert pytest.approx(err) == 15 * [1.0] + + assert minpack.library.lmder1(fcn, x, fvec, fjac, tol) == 1 + + assert pytest.approx(x, abs=100 * tol) == [0.8241058e-1, 0.1133037e1, 0.2343695e1] + + +def test_lmder1_exception(): + + class DummyException(Exception): ... + + def fcn(x, fvec, fjac, jacobian: bool) -> None: + raise DummyException() + + x = np.array([1.0, 1.0, 1.0]) + fvec = np.zeros(15, dtype=np.float64) + fjac = np.zeros((3, 15), dtype=np.float64) + tol = sqrt(np.finfo(np.float64).eps) + + with pytest.raises(DummyException): + minpack.library.lmder1(fcn, x, fvec, fjac, tol) + + +def test_lmdif1(): + y = np.array( + [ + 1.4e-1, + 1.8e-1, + 2.2e-1, + 2.5e-1, + 2.9e-1, + 3.2e-1, + 3.5e-1, + 3.9e-1, + 3.7e-1, + 5.8e-1, + 7.3e-1, + 9.6e-1, + 1.34e0, + 2.1e0, + 4.39e0, + ] + ) + + def fcn(x, fvec) -> None: + for i in range(fvec.size): + tmp1, tmp2 = i + 1, 16 - i - 1 + tmp3 = tmp2 if i >= 8 else tmp1 + fvec[i] = y[i] - (x[0] + tmp1 / (x[1] * tmp2 + x[2] * tmp3)) + + x = np.array([1.0, 1.0, 1.0]) + fvec = np.zeros(15, dtype=np.float64) + fjac = np.zeros((3, 15), dtype=np.float64) + tol = sqrt(np.finfo(np.float64).eps) + + assert minpack.library.lmdif1(fcn, x, fvec, tol) == 1 + + assert pytest.approx(x, abs=100 * tol) == [0.8241058e-1, 0.1133037e1, 0.2343695e1] + + +def test_lmdif1_exception(): + + class DummyException(Exception): ... + + def fcn(x, fvec) -> None: + raise DummyException() + + x = np.array([1.0, 1.0, 1.0]) + fvec = np.zeros(15, dtype=np.float64) + fjac = np.zeros((3, 15), dtype=np.float64) + tol = sqrt(np.finfo(np.float64).eps) + + with pytest.raises(DummyException): + minpack.library.lmdif1(fcn, x, fvec, tol) diff --git a/python/setup.cfg b/python/setup.cfg new file mode 100644 index 0000000..e47d585 --- /dev/null +++ b/python/setup.cfg @@ -0,0 +1,38 @@ +[metadata] +name = minpack +version = 2.0.0 +desciption = Python bindings for MINPACK +long_desciption = file: README.rst +long_description_content_type = text/x-rst +url = https://github.com/fortran-lang/minpack +license = MIT +license_files = + ../LICENSE.txt +classifiers = + Development Status :: 5 - Production + Intended Audience :: Science/Research + Programming Language :: Fortran + Programming Language :: Python :: 3 :: Only + Programming Language :: Python :: 3 + Programming Language :: Python :: 3.6 + Programming Language :: Python :: 3.7 + Programming Language :: Python :: 3.8 + Programming Language :: Python :: 3.9 + Programming Language :: Python :: 3.10 + +[options] +packages = find: +install_requires = + cffi + numpy +tests_require = + pytest + pytest-cov +python_requires = >=3.6 + +[coverage:run] +omit = + */test_*.py + +[aliases] +test = pytest diff --git a/python/setup.py b/python/setup.py new file mode 100644 index 0000000..f48c3eb --- /dev/null +++ b/python/setup.py @@ -0,0 +1,7 @@ +from setuptools import setup + +setup( + cffi_modules=["ffi-builder.py:ffibuilder"], + package_data={"minpack": ["_libminpack*.so"]}, +) + diff --git a/test/api/tester.c b/test/api/tester.c index 03fb536..d609490 100644 --- a/test/api/tester.c +++ b/test/api/tester.c @@ -210,18 +210,18 @@ test_lmder1 (void) minpack_chkder(m, n, x, fvec, fjac, m, xp, fvecp, 1, err); info = 1; - trial_lmder_fcn(m, n, x, fvec, fjac, m, &info, y); + trial_lmder_fcn(m, n, x, fvec, fjac, m, &info, (void*)y); info = 2; - trial_lmder_fcn(m, n, x, fvec, fjac, m, &info, y); + trial_lmder_fcn(m, n, x, fvec, fjac, m, &info, (void*)y); info = 1; - trial_lmder_fcn(m, n, xp, fvecp, fjac, m, &info, y); + trial_lmder_fcn(m, n, xp, fvecp, fjac, m, &info, (void*)y); minpack_chkder(m, n, x, fvec, fjac, m, xp, fvecp, 2, err); for (int i = 0; i < 15; i++) { if (!check(err[i], 1.0, tol, "Unexpected derivatives")) return 1; } - minpack_lmder1(trial_lmder_fcn, m, n, x, fvec, fjac, m, tol, &info, ipvt, wa, 30, y); + minpack_lmder1(trial_lmder_fcn, m, n, x, fvec, fjac, m, tol, &info, ipvt, wa, 30, (void*)y); if (!check(info, 1, "Unexpected info value")) return 1; if (!check(x[0], 0.8241058e-1, 100*tol, "Unexpected x[0]")) return 1; if (!check(x[1], 0.1133037e+1, 100*tol, "Unexpected x[1]")) return 1; @@ -247,11 +247,11 @@ test_lmder (void) minpack_chkder(m, n, x, fvec, fjac, m, xp, fvecp, 1, err); info = 1; - trial_lmder_fcn(m, n, x, fvec, fjac, m, &info, y); + trial_lmder_fcn(m, n, x, fvec, fjac, m, &info, (void*)y); info = 2; - trial_lmder_fcn(m, n, x, fvec, fjac, m, &info, y); + trial_lmder_fcn(m, n, x, fvec, fjac, m, &info, (void*)y); info = 1; - trial_lmder_fcn(m, n, xp, fvecp, fjac, m, &info, y); + trial_lmder_fcn(m, n, xp, fvecp, fjac, m, &info, (void*)y); minpack_chkder(m, n, x, fvec, fjac, m, xp, fvecp, 2, err); for (int i = 0; i < 15; i++) { @@ -259,7 +259,7 @@ test_lmder (void) } minpack_lmder(trial_lmder_fcn, m, n, x, fvec, fjac, m, tol, tol, 0.0, 2000, diag, 1, - 100.0, 0, &info, &nfev, &njev, ipvt, qtf, wa1, wa2, wa3, wa4, y); + 100.0, 0, &info, &nfev, &njev, ipvt, qtf, wa1, wa2, wa3, wa4, (void*)y); if (!check(info, 1, "Unexpected info value")) return 1; if (!check(x[0], 0.8241058e-1, 100*tol, "Unexpected x[0]")) return 1; if (!check(x[1], 0.1133037e+1, 100*tol, "Unexpected x[1]")) return 1; @@ -298,7 +298,7 @@ test_lmdif1 (void) int lwa = m*n + 5*n + m; double wa[lwa]; - minpack_lmdif1(trial_lmdif_fcn, 15, 3, x, fvec, tol, &info, ipvt, wa, lwa, y); + minpack_lmdif1(trial_lmdif_fcn, 15, 3, x, fvec, tol, &info, ipvt, wa, lwa, (void*)y); if (!check(info, 1, "Unexpected info value")) return 1; if (!check(x[0], 0.8241058e-1, 100*tol, "Unexpected x[0]")) return 1; if (!check(x[1], 0.1133037e+1, 100*tol, "Unexpected x[1]")) return 1; @@ -321,7 +321,7 @@ test_lmdif (void) double fjac[m*n], diag[n], qtf[n], wa1[n], wa2[n], wa3[n], wa4[m]; minpack_lmdif(trial_lmdif_fcn, 15, 3, x, fvec, tol, tol, 0.0, 2000, 0.0, diag, 1, - 100.0, 0, &info, &nfev, fjac, 15, ipvt, qtf, wa1, wa2, wa3, wa4, y); + 100.0, 0, &info, &nfev, fjac, 15, ipvt, qtf, wa1, wa2, wa3, wa4, (void*)y); if (!check(info, 1, "Unexpected info value")) return 1; if (!check(x[0], 0.8241058e-1, 100*tol, "Unexpected x[0]")) return 1; if (!check(x[1], 0.1133037e+1, 100*tol, "Unexpected x[1]")) return 1; diff --git a/test/meson.build b/test/meson.build index 6279d4b..1f744e6 100644 --- a/test/meson.build +++ b/test/meson.build @@ -19,7 +19,7 @@ foreach t : tests ) endforeach -if add_languages('c', required: false, native: false) +if has_cc test( 'c-api', executable(