jonathf / numpoly

Numpy compatible polynomial representation
https://numpoly.readthedocs.io
BSD 2-Clause "Simplified" License
12 stars 6 forks source link

Addition of array of polynomials with scalar raises ValueError when comparing exponents. #101

Open JulVandenBroeck opened 1 year ago

JulVandenBroeck commented 1 year ago

Hello jonathf!

I have been using your modules in my research for a while now and I love them. I have stumbled upon a problem I can't manage to solve, however.

I have an array of polynomials (with 2 elements in the example) and want to add 1 to all of them. Running the following code:

import numpoly
qs = numpoly.variable(2)
qs+1

Raises the following error:

Traceback (most recent call last):
  File "/home/quest/juvdnbro/repos/gpc/src/schrodinger1d/test.py", line 3, in <module>
    qs+1
  File "/home/quest/juvdnbro/.local/lib/python3.10/site-packages/numpoly/baseclass.py", line 242, in __array_ufunc__
    return numpoly.UFUNC_COLLECTION[ufunc](*inputs, **kwargs)
  File "/home/quest/juvdnbro/.local/lib/python3.10/site-packages/numpoly/array_function/add.py", line 60, in add
    return simple_dispatch(
  File "/home/quest/juvdnbro/.local/lib/python3.10/site-packages/numpoly/dispatch.py", line 84, in simple_dispatch
    inputs = numpoly.align_polynomials(*inputs)
  File "/home/quest/juvdnbro/.local/lib/python3.10/site-packages/numpoly/align.py", line 47, in align_polynomials
    polys = align_exponents(*polys)
  File "/home/quest/juvdnbro/.local/lib/python3.10/site-packages/numpoly/align.py", line 211, in align_exponents
    if numpy.all(poly.exponents == global_exponents):
ValueError: operands could not be broadcast together with shapes (2,2) (3,2)

It seems that the scalar 1 gets converted into polynomial([1,1]), but the exponents don't have the correct shape. Is there a better way to do this, or would it be possible to allow for this kind of arithmetic?

Thanks in advance! Jul

JulVandenBroeck commented 1 year ago

I have found that the exponents of qs are [array([1, 0]), array([0, 1])], and if I do the following hack (adding 1 and subtracting again):

qs_hack = numpoly.polynomial([q0+1,q1+1]) - 1

the exponents become [array([0, 0]), array([1, 0]), array([0, 1])]. This form allows for the arithmetic that I require. Is there a more elegant way for this? 😆

JulVandenBroeck commented 1 year ago

Another update: with this hack, division by scalars doesn't work. qs_hack/2 raises the same exception.

Simply commenting out the following lines at the end of align.py solves my problem:

# if numpy.all(poly.exponents == global_exponents):
#     continue
jonathf commented 1 year ago

I'm not reproducing the error on my end. Which version of numpoly are you on?

Also, can you compare version v1.2.9 and v1.2.10 and see if that makes a difference?

JulVandenBroeck commented 1 year ago

Thanks for the reply. I was using v1.2.9 with Python 3.10.6. I tried both v1.2.9 and v1.2.10 now and for me the error occurs for both versions... Any ideas?

JulVandenBroeck commented 1 year ago

Also, I'm using Numpy v1.25.0.

jonathf commented 1 year ago

Updating numpy to v1.25 did the trick; I can reproduce the issue. I'll look into the issue when I have a moment to spare.

jonathf commented 1 year ago

Numpy 1.25 had a few incompatabilities., but the issue should now be resolved in numpoly v1.2.11. Let me know if there are any issues.

JulVandenBroeck commented 1 year ago

It works perfectly now. Great work, I really appreciate it!