A permissively licensed modern implementation of Gauss-Jacobi quadrature which returns the weights and nodes over the standard interval [-1, 1].
The most automated approach is to use the conda environment and fpm build:
micromamba create -f environment.yml
micromamba activate gaussjacquad
fpm buildAn analytic result can be obtained from the scripts folder.
python scripts/sympy_gauss_jac.py --npts <npoints> --alpha <alpha> --beta <beta> --n_dig <precision>
fpm testRunning the implemented recursion based Gauss-Jacobi can be done via:
fpm run gjp_quad_rec -- <npoints> <alpha> <beta>
fpm run gjp_quad -- <npoints> <alpha> <beta> <method>Currently the only supported methods are “rec” and “gw” (Golub Welsch).
A meson build backend is also present, which makes it easier to incorporate as subprojects of projects other than those supported by fpm.
meson setup bbdir
meson compile -C bbdir
./bbdir/gjp_quad <npoints> <alpha> <beta> <method>We support and encourage users to generate single file versions of the algorithms here to include in their code bases. This can be done with the scripts/add_headers.py script:
# Strips comments by default
python scripts/export_single.py --modulename "gjp_algo665"
python scripts/export_single.py --modulename "gjp_gw" --keep-commentsThese can be dropped into any code base or compiled as is into a shared library.
gfortran dist/gjp_algo665_single.f90 --shared -fPIC -o libgjp_algo665.soWe provide a header only interface, which bypasses passing strings between
Fortran and C. In order to do this, there is some duplication logic in
GaussJacobiQuad and GaussJacobiQuadCInterp.
There is also a CLI interface to the C bound interface, c_cli_gjpq. However,
this will not be compiled by fpm.
meson setup bbdir
meson compile -C bbdir
./bbdir/c_cli_gjpq <npoints> <alpha> <beta> <method>The C CLI might be more pleasant in that decimals do not need to be provided
explicitly for alpha and beta.
The ISO_C_BINDING variant of the code is also used for a python
interface generated with f2py. It is easiest to use with the new meson back-end:
cd interfaces/PyInterface
f2py -c --backend meson gjquadpy.pyf \
--dep lapack \
../../src/GaussJacobiQuadCCompat.f90 \
../../src/GaussJacobiQuad.f90 \
../../src/gjp_constants.f90 \
../../src/gjp_types.f90 \
../../src/gjp_rec.f90 \
../../src/gjp_common.f90 \
../../src/gjp_lapack.f90 \
../../src/gjp_gw.f90 \
../../src/gjp_algo665.f90Once compiled, the gjpquad_cli.py script can be used to run the code:
python gjquad_cli.py --npts 5 --alpha 2 --beta 3
Root: -6.90457750126761027E-01 Weight: 2.74101780663370022E-02
Root: -3.26519931349000647E-01 Weight: 2.12917860603648035E-01
Root: 8.23378495520349085E-02 Weight: 4.39084379443950568E-01
Root: 4.75178870612831761E-01 Weight: 3.22206565472217876E-01
Root: 7.92794294644228348E-01 Weight: 6.50476830805121059E-02- The
recmethod fails for high values ofbetaso thegwmethod
should be used in such situations.
algo665is an in-place variant ofgwand is much faster when many points are needed
A very preliminary set can be run once all the interfaces have been compiled:
# From $GITROOT
mv interfaces/PyInterface/gjquadpy*.so .
hyperfine --warmup 10 --export-markdown gjp_benchmarks.md \
'bbdir/gjp_quad 5 2. 3. "gw"' \
'bbdir/c_cli_gjpq 5 2 3 gw' \
'PYTHONPATH=$(pwd) python interfaces/PyInterface/gjquad_cli.py --npts 5 --alpha 2 --beta 3' \
'python scripts/scipy_gauss_jac.py --npts 5 --alpha 2 --beta 3' \
'python scripts/sympy_gauss_jac.py --npts 5 --alpha 2 --beta 3 --n_dig 15'Which gives:
ummary
bbdir/c_cli_gjpq 5 2 3 gw ran
1.03 ± 0.45 times faster than bbdir/gjp_quad 5 2. 3. "gw"
39.08 ± 11.48 times faster than PYTHONPATH=$(pwd) python interfaces/PyInterface/gjquad_cli.py --npts 5 --alpha 2 --beta 3
55.71 ± 16.51 times faster than python scripts/scipy_gauss_jac.py --npts 5 --alpha 2 --beta 3
86.73 ± 24.86 times faster than python scripts/sympy_gauss_jac.py --npts 5 --alpha 2 --beta 3 --n_dig 15
hyperfine --warmup 10 --export-markdown gjp_benchmarks.md 32.65s user 61.31s system 469% cpu 20.004 totalOr in other words:
| Command | Mean [ms] | Min [ms] | Max [ms] | Relative |
|---|---|---|---|---|
gjp_quad 5 2. 3. "gw" | 2.6 ± 0.8 | 1.8 | 7.8 | 1.03 ± 0.45 |
c_cli_gjpq 5 2 3 gw | 2.5 ± 0.7 | 1.8 | 8.4 | 1.00 |
gjquad_cli.py --npts 5 --alpha 2 --beta 3 | 97.5 ± 6.3 | 95.4 | 130.8 | 39.08 ± 11.48 |
scipy_gauss_jac.py --npts 5 --alpha 2 --beta 3 | 139.0 ± 10.6 | 132.4 | 173.4 | 55.71 ± 16.51 |
sympy_gauss_jac.py --npts 5 --alpha 2 --beta 3 --n_dig 15 | 216.4 ± 1.7 | 214.7 | 219.7 | 86.73 ± 24.86 |
Which suggests what one might suspect, that there is a large overhead in calling
python , and that the C and Fortran variants are almost exactly as fast as
each other. However, the f2py variant is still way faster than the existing
python implementations.
hyperfine --warmup 10 --export-markdown gjp_benchmarks.md \
'PYTHONPATH=$(pwd) python interfaces/PyInterface/gjquad_cli.py --npts 5 --alpha 2 --beta 3' \
'python scripts/scipy_gauss_jac.py --npts 5 --alpha 2 --beta 3'
Summary
PYTHONPATH=$(pwd) python interfaces/PyInterface/gjquad_cli.py --npts 5 --alpha 2 --beta 3 ran
1.38 ± 0.11 times faster than python scripts/scipy_gauss_jac.py --npts 5 --alpha 2 --beta 3A pre-commit job is setup on CI to enforce consistent styles, so it is best to
set it up locally as well (using pipx for isolation):
# Run before commiting
pipx run pre-commit run --all-files
# Or install the git hook to enforce this
pipx run pre-commit installWhen the headers in the sources need to be updated modify add_headers.py and run:
python scripts/add_headers.py --dirs src/ interfaces/ --ftypes "f90,f77" --cchar '!'
python scripts/add_headers.py --dirs interfaces --ftypes "c,h" --cchar '//'
python scripts/add_headers.py --dirs interfaces scripts --ftypes "py" --cchar '#'Remember to do this before exporting the code into other projects (e.g. featom).
MIT.
If you use this library (including the interfaces) please remember to cite it as:
Rohit Goswami. (2023). HaoZeke/GaussJacobiQuad: GaussJacobiQuad I (v0.1.0). Zenodo. https://doi.org/10.5281/ZENODO.8285112
Or use the bibtex entry:
@software{rgGaussQuad23,
author = {Rohit Goswami},
title = {HaoZeke/GaussJacobiQuad: GaussJacobiQuad I},
month = aug,
year = 2023,
publisher = {Zenodo},
version = {v0.1.0},
doi = {10.5281/zenodo.8285112},
url = {https://doi.org/10.5281/zenodo.8285112}
}An ArXiv –> JOSS paper is in the works.