search

jonathf / chaospy

Chaospy - Toolbox for performing uncertainty quantification.
https://chaospy.readthedocs.io/
MIT License
439 stars 87 forks source link

High-Order 1-Dimensional Gaussian Quadrature Bug #420

Open ehsansaleh opened 5 months ago

ehsansaleh commented 5 months ago

Describe the bug 1-dimensional Gaussian quadrature with orders more than 118 throws out a numpoly error.

To Reproduce

>>> import chaospy

>>> chaospy.__version__
'4.3.13'

>>> cpy_dist = chaospy.Iid(chaospy.Uniform(0, 1), 1)
>>> ex_quadx_np, ex_quadw_np = chaospy.generate_quadrature(
...     order=128, dist=cpy_dist, rule='gaussian', 
...     recurrence_algorithm='stieltjes', sparse=False)

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "./venv/lib/python3.8/site-packages/chaospy/quadrature/frontend.py", line 171, in generate_quadrature
    abscissas, weights = _generate_quadrature(
  File "./venv/lib/python3.8/site-packages/chaospy/quadrature/frontend.py", line 287, in _generate_quadrature
    abscissas, weights = quad_function(order, dist, **parameters)
  File "./venv/lib/python3.8/site-packages/chaospy/quadrature/gaussian.py", line 81, in gaussian
    coefficients = chaospy.construct_recurrence_coefficients(
  File "./venv/lib/python3.8/site-packages/chaospy/recurrence/frontend.py", line 117, in construct_recurrence_coefficients
    coeffs, _, _ = stieltjes(order, dist, rule=rule, tolerance=tolerance)
  File "./venv/lib/python3.8/site-packages/chaospy/recurrence/stieltjes.py", line 68, in stieltjes
    return analytical_stieltjes(order=order, dist=dist)
  File "./venv/lib/python3.8/site-packages/chaospy/recurrence/stieltjes.py", line 180, in analytical_stieltjes
    orth = numpoly.polynomial(orth[1:]).T
  File "./venv/lib/python3.8/site-packages/numpoly/construct/polynomial.py", line 137, in polynomial
    poly = compose_polynomial_array(
  File "./venv/lib/python3.8/site-packages/numpoly/construct/compose.py", line 41, in compose_polynomial_array
    [
  File "./venv/lib/python3.8/site-packages/numpoly/construct/compose.py", line 42, in <listcomp>
    numpoly.expand_dims(numpoly.aspolynomial(array, dtype=dtype), axis=0)
  File "./venv/lib/python3.8/site-packages/numpoly/array_function/expand_dims.py", line 38, in expand_dims
    return numpoly.polynomial(out, names=a.indeterminants)
  File "./venv/lib/python3.8/site-packages/numpoly/construct/polynomial.py", line 104, in polynomial
    poly = numpoly.ndpoly.from_attributes(
  File "./venv/lib/python3.8/site-packages/numpoly/baseclass.py", line 387, in from_attributes
    return numpoly.polynomial_from_attributes(
  File "./venv/lib/python3.8/site-packages/numpoly/construct/from_attributes.py", line 72, in polynomial_from_attributes
    exponents, coefficients, names = clean.postprocess_attributes(
  File "./venv/lib/python3.8/site-packages/numpoly/construct/clean.py", line 106, in postprocess_attributes
    raise PolynomialConstructionError(
numpoly.construct.clean.PolynomialConstructionError: expected len(exponents) == len(coefficients_); found 119 != 120
>>>

Expected behavior It should not raise any errors. This is concerning since smaller orders may be suffering from a related bug.

Desktop (please complete the following information):

Additional context None

ehsansaleh commented 5 months ago

Here is my pip freeze output if it helps:

absl-py==1.4.0
anyio==4.3.0
argon2-cffi==23.1.0
argon2-cffi-bindings==21.2.0
arrow==1.3.0
asttokens==2.2.1
astunparse==1.6.3
async-lru==2.0.4
attrs==22.2.0
Babel==2.14.0
backcall==0.2.0
beautifulsoup4==4.12.3
bleach==6.1.0
blosc2==2.0.0
bokeh==3.0.3
cachetools==5.3.0
certifi==2022.12.7
cffi==1.16.0
chaospy==4.3.13
charset-normalizer==3.0.1
comm==0.1.2
contourpy==1.0.7
cycler==0.11.0
Cython==0.29.33
debugpy==1.6.5
decorator==5.1.1
deepdiff==6.2.3
defusedxml==0.7.1
entrypoints==0.4
exceptiongroup==1.2.0
executing==1.2.0
fastjsonschema==2.16.2
filelock==3.13.4
flatbuffers==23.1.21
fonttools==4.38.0
fqdn==1.5.1
fsspec==2024.3.1
gast==0.4.0
google-auth==2.16.1
google-auth-oauthlib==0.4.6
google-pasta==0.2.0
grpcio==1.51.1
h11==0.14.0
h5py==3.7.0
httpcore==1.0.5
httpx==0.27.0
idna==3.4
importlib-metadata==6.0.0
importlib-resources==5.10.2
ipykernel==6.20.2
ipython==8.8.0
ipython-genutils==0.2.0
isoduration==20.11.0
jedi==0.18.2
Jinja2==3.1.2
json5==0.9.25
jsonpointer==2.4
jsonschema==4.21.1
jsonschema-specifications==2023.12.1
jupyter-contrib-core==0.4.2
jupyter-events==0.10.0
jupyter-lsp==2.2.5
jupyter-nbextensions-configurator==0.6.3
jupyter_client==7.4.9
jupyter_core==5.1.3
jupyter_server==2.14.0
jupyter_server_terminals==0.5.3
jupyterlab==4.1.6
jupyterlab_pygments==0.3.0
jupyterlab_server==2.26.0
jupytext==1.14.4
keras==2.11.0
kiwisolver==1.4.4
libclang==15.0.6.1
Markdown==3.4.1
markdown-it-py==2.1.0
MarkupSafe==2.1.2
matplotlib==3.6.3
matplotlib-inline==0.1.6
mdit-py-plugins==0.3.3
mdurl==0.1.2
mistune==3.0.2
mpmath==1.3.0
msgpack==1.0.4
nbclient==0.10.0
nbconvert==7.16.3
nbformat==5.7.3
nest-asyncio==1.5.6
networkx==3.1
notebook==6.4.12
notebook_shim==0.2.4
numexpr==2.8.4
numpoly==1.2.11
numpy==1.24.1
nvidia-cublas-cu12==12.1.3.1
nvidia-cuda-cupti-cu12==12.1.105
nvidia-cuda-nvrtc-cu12==12.1.105
nvidia-cuda-runtime-cu12==12.1.105
nvidia-cudnn-cu12==8.9.2.26
nvidia-cufft-cu12==11.0.2.54
nvidia-curand-cu12==10.3.2.106
nvidia-cusolver-cu12==11.4.5.107
nvidia-cusparse-cu12==12.1.0.106
nvidia-nccl-cu12==2.19.3
nvidia-nvjitlink-cu12==12.4.127
nvidia-nvtx-cu12==12.1.105
oauthlib==3.2.2
opt-einsum==3.3.0
ordered-set==4.1.0
orjson==3.8.5
overrides==7.7.0
packaging==23.0
pandas==1.5.2
pandocfilters==1.5.1
parso==0.8.3
pexpect==4.8.0
pickle5==0.0.11
pickleshare==0.7.5
Pillow==9.4.0
pkgutil_resolve_name==1.3.10
platformdirs==2.6.2
prometheus_client==0.20.0
prompt-toolkit==3.0.36
protobuf==3.19.6
psutil==5.9.4
ptyprocess==0.7.0
pure-eval==0.2.2
py-cpuinfo==9.0.0
pyasn1==0.4.8
pyasn1-modules==0.2.8
pycparser==2.22
Pygments==2.14.0
pyinstrument==4.4.0
pyparsing==3.0.9
pyrsistent==0.19.3
python-dateutil==2.8.2
python-json-logger==2.0.7
pytz==2022.7.1
PyYAML==6.0
pyzmq==25.1.2
referencing==0.34.0
requests==2.31.0
requests-oauthlib==1.3.1
rfc3339-validator==0.1.4
rfc3986-validator==0.1.1
rpds-py==0.18.0
rsa==4.9
ruamel.yaml==0.17.21
ruamel.yaml.clib==0.2.7
scipy==1.10.0
seaborn==0.12.2
Send2Trash==1.8.3
six==1.16.0
sniffio==1.3.1
soupsieve==2.5
stack-data==0.6.2
sympy==1.12
tables==3.8.0
tensorboard==2.11.2
tensorboard-data-server==0.6.1
tensorboard-plugin-wit==1.8.1
tensorboardX==2.5.1
tensorflow==2.11.0
tensorflow-estimator==2.11.0
tensorflow-io-gcs-filesystem==0.30.0
termcolor==2.2.0
terminado==0.18.1
tinycss2==1.2.1
toml==0.10.2
tomli==2.0.1
torch==2.2.2
torchvision==0.17.2
tornado==6.2
traitlets==5.9.0
triton==2.2.0
types-python-dateutil==2.9.0.20240316
typing_extensions==4.11.0
uri-template==1.3.0
urllib3==1.26.14
wcwidth==0.2.6
webcolors==1.13
webencodings==0.5.1
websocket-client==1.7.0
Werkzeug==2.2.3
wrapt==1.14.1
xyzservices==2023.2.0
zipp==3.11.0
jonathf commented 4 months ago

You reached the limit of how big one exponent can be. I've updated in version 4.3.15. Let me know if that solves the issue for you.

EDIT: You also need to update numpoly to version 1.2.12.