About

A permissively licensed modern implementation of Gauss-Jacobi quadrature which returns the weights and nodes over the standard interval [-1, 1].

Usage

The most automated approach is to use the conda environment and fpm build:

micromamba create -f environment.yml
micromamba activate gaussjacquad
fpm build

An 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 test

Running 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).

Meson Support

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>

External Fortran libaries

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-comments

These 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.so

Interfaces

C/C++ header interface

We 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.

f2py generated interface

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.f90

Once 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

Notes

The rec method fails for high values of beta so the gw method

should be used in such situations.

algo665 is an in-place variant of gw and is much faster when many points are needed

Benchmarks

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 total

Or 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 3

Development

Developing locally

A 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 install

Updating licenses

When 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).

License

MIT.

Citation

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.