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

Code

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

Code

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>

Code

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>

Code

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

Code

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

Code

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>

Code

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

Code

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

Code

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'

Code

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

Code

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

Code

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

Code

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 '#'

Code

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

Code

An ArXiv –> JOSS paper is in the works.