jonathf / numpoly

Numpy compatible polynomial representation

Hello there!

I'm Jonathan, and I am an all-around nerd.

This is most observable online through my dotfiles. There you can see I like to tinker around with awesome wm, qutebrowser, and neovim, among other things. I also have an above average interest in gin, grabbing about 95 % of my liquor cabinet space. I also would say I have hobbies, but after children started to occupy my living room, my time to tinker with my various projects has been reduced a bit.

I have a PhD in applied mathematics and probability theory where I learning a few things about polynomial chaos expansions. I wrote a Python library about it called chaospy, which some people seem to like. I still maintain it and add some features there now and again. I also maintain pyvroom a library for solving solving various vehicle routing problems.

I love doing research, and if it weren't for the tedious paper and funding application writing, I'd likely still be doing that. Instead I work for Expert Analytics, where I mostly work as a software developer and data scientist, and have a lot fun.

Here are some stuff I think are neat:

numpoly's Issues

Diagonalize parametric matrix

Hello,

is it possible to get a parametric eigendecomposition of a matrix whose elements are a mix of known numerical values and parametric variables, expressed as ndpoly entries?

Thank you.

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

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

Min/max/argmin/argmax support

• For polynomials, first sort lexicographically by largest exponent: x**3 > x*y > y
• For equal order, sort by largest lead coefficient: 4*x**2 > 2*x*y > -8*y**2

This alow for being backward compatible with numpy sorting methods as constant is a lead coefficient.

AttributeError with division

Hi,

I'm a brand new user on Github and therefore sorry if it is not how I'm expected to use this topic.

Your package seems very useful and I thank you for implementing all of that. However, I'm finding an error when dividing a polynomial by another. Using "/", "%" and "divmod" all end up in an error :
AttributeError: module 'numpy' has no attribute 'broadcast_shapes'

Do you have any idea what in the division functions could end up in such an error ? I tried with elements of the 'numpoly.baseclass.ndpoly' class of different lengths, and none worked. Note that the same error appears when using the symbol "*".

Thank you,
Archer

Make proxy array function

The features suggested in #54 depend on sorting the polynomial. Instead of sorting in place, create a function that substitute polynomial for integers that whose order correspond to the sorting in question.

Easies way to solve #54 for now.

function for gathering all lead exponents

Same as #53, but returns the actual exponent instead of the coefficient.

six is not in install_requires but numpoly/baseclass.py imports six

Problem with flatiter

I believe the flatiter might not be handled correctly.

import numpoly
X = numpoly.variable(6).reshape((3,2))
print('using flatten:')
for p in X.flatten():
    print(p)
print('using flat:')
for p in X.flat:
    print(p)

output:

using flatten:
q0
q1
q2
q3
q4
q5
using flat:
(1, 0, 0, 0, 0, 0)
(0, 1, 0, 0, 0, 0)
(0, 0, 1, 0, 0, 0)
(0, 0, 0, 1, 0, 0)
(0, 0, 0, 0, 1, 0)
(0, 0, 0, 0, 0, 1)

``power`` implementation

Current implementation is a little hacky: repeat multiplying against itself n
times. For-loop over n, if array.

• Create a method that does not rely on repeated calling multiply, and instead
  allocates a single chunk of memory and fill inn results there.
• Somehow avoid using element-by-element for-loop over exponents.

Add support save/load to disk

Pickle should work out of the box. Preferably using the numpy.save
interface.

test failed

test/test_align.py . [ 0%]
test/test_array_function.py ............................................ [ 35%]
....................................................................... [ 90%]
test/test_baseclass.py F [ 91%]
test/test_construct.py .. [ 92%]
test/test_poly_function.py ......... [100%]
=================================== FAILURES ===================================
_________________________________ test_scalars _________________________________
def test_scalars():
assert XY.shape == (2,)
assert XY.size == 2
assert X.shape == ()
assert X.size == 1
assert EMPTY.shape in [(), (0,)] # different behavior in py2/3
assert EMPTY.size == 0

    assert numpy.all(numpy.array(XY.coefficients) == [[1, 0], [0, 1]])
    assert X.coefficients == [1]
    assert EMPTY.coefficients == []

    assert numpy.all(XY.exponents == [[1, 0], [0, 1]])
    assert XY.exponents.shape == (2, 2)
    assert X.exponents == 1
    assert X.exponents.shape == (1, 1)
    assert numpy.all(EMPTY.exponents == 0)
    assert EMPTY.exponents.shape == (1, 1)

    assert numpy.all(XY.indeterminants == XY)
    assert X.indeterminants == X

    assert numpy.all(XY.values == numpy.array(
        [(1, 0), (0, 1)], dtype=[("<;", "<i8",), (";<", "<i8",)]))
    assert X.values == numpy.array((1,), dtype=[("<", "<i8",)])

  assert EMPTY.values.dtype == numpy.dtype([(";", "<i8")])

E AssertionError: assert dtype([(';', '<i4')]) == dtype([(';', '<i8')])
E + where dtype([(';', '<i4')]) = array([], dtype=[(';', '<i4')]).dtype
E + where array([], dtype=[(';', '<i4')]) = EMPTY.values
E + and dtype([(';', '<i8')]) = <class 'numpy.dtype'>([(';', '<i8')])
E + where <class 'numpy.dtype'> = numpy.dtype
test/test_baseclass.py:33: AssertionError
=============================== warnings summary ===============================
test/test_array_function.py::test_numpy_isfinite[func_interface0]
test/test_array_function.py::test_numpy_isfinite[func_interface1]
/builddir/build/BUILD/numpoly-0.1.16/test/test_array_function.py:236: RuntimeWarning: invalid value encountered in log
poly = numpoly.polynomial([numpy.log(-1.), X, numpy.log(0)])
test/test_array_function.py::test_numpy_isfinite[func_interface0]
test/test_array_function.py::test_numpy_isfinite[func_interface1]
/builddir/build/BUILD/numpoly-0.1.16/test/test_array_function.py:236: RuntimeWarning: divide by zero encountered in log
poly = numpoly.polynomial([numpy.log(-1.), X, numpy.log(0)])
test/test_array_function.py::test_numpy_nonzero[interface0]
test/test_array_function.py::test_numpy_nonzero[interface1]
/usr/lib/python3.8/site-packages/numpy/core/fromnumeric.py:61: DeprecationWarning: Calling nonzero on 0d arrays is deprecated, as it behaves surprisingly. Use atleast_1d(cond).nonzero() if the old behavior was intended. If the context of this warning is of the form arr[nonzero(cond)], just use arr[cond].
return bound(*args, **kwds)
test/test_array_function.py::test_numpy_nonzero[interface2]
/builddir/build/BUILD/numpoly-0.1.16/conftest.py:17: DeprecationWarning: Calling nonzero on 0d arrays is deprecated, as it behaves surprisingly. Use atleast_1d(arr).nonzero() if the old behavior was intended.
return func(*args, **kwargs)
test/test_poly_function.py::test_numpoly_cross_truncate
/builddir/build/BUILD/numpoly-0.1.16/numpoly/poly_function/monomial/cross_truncation.py:49: RuntimeWarning: divide by zero encountered in true_divide
out = numpy.sum((indices/bound)norm, axis=-1)(1./norm) <= 1
test/test_poly_function.py::test_numpoly_cross_truncate
/builddir/build/BUILD/numpoly-0.1.16/numpoly/poly_function/monomial/cross_truncation.py:49: RuntimeWarning: invalid value encountered in true_divide
out = numpy.sum((indices/bound)norm, axis=-1)(1./norm) <= 1
test/test_poly_function.py::test_numpoly_cross_truncate
/builddir/build/BUILD/numpoly-0.1.16/numpoly/poly_function/monomial/cross_truncation.py:49: RuntimeWarning: invalid value encountered in less_equal
out = numpy.sum((indices/bound)norm, axis=-1)(1./norm) <= 1
test/test_poly_function.py::test_numpoly_cross_truncate
/builddir/build/BUILD/numpoly-0.1.16/numpoly/poly_function/monomial/cross_truncation.py:47: RuntimeWarning: divide by zero encountered in true_divide
out = numpy.max(indices/bound, axis=-1) <= 1
test/test_poly_function.py::test_numpoly_cross_truncate
/builddir/build/BUILD/numpoly-0.1.16/numpoly/poly_function/monomial/cross_truncation.py:47: RuntimeWarning: invalid value encountered in true_divide
out = numpy.max(indices/bound, axis=-1) <= 1
test/test_poly_function.py::test_numpoly_cross_truncate
/builddir/build/BUILD/numpoly-0.1.16/numpoly/poly_function/monomial/cross_truncation.py:47: RuntimeWarning: invalid value encountered in less_equal
out = numpy.max(indices/bound, axis=-1) <= 1
-- Docs: https://docs.pytest.org/en/latest/warnings.html
============== 1 failed, 127 passed, 13 warnings in 15.83 seconds ==============

Please update the version on the PyPI website

here

Preserving indeterminants with numpoly.multiply etc

Is there a way to preserve the indeterminants of a set of polynomials when one applies multiply, sum, etc? For example, I create a set of 3-dimensional polynomials which may contain indeterminants q0, q1, q2 but sometimes only q0 (like when you create a basis of graded lexicographic monomials in d-dimensions). After applying align_indeterminants() I now have a set where all polynomials return ('q0, 'q1', 'q2') as "names" even though a polynomial may only be 'q0^2' etc. However, when I use numpoly.multiply(q0^2, q1) I get a new polynomial that has as names ('q0', 'q1') only. I would like to preserve the structure I had before where the new polynomial still returns ('q0, 'q1', 'q2') as "names" even though indeterminate 'q2' does not appear in the polynomial. This is useful because I have a dictionary that has as keys the exponents in the form [., ., .].

Element-in support: ``x in y`` and ``{numpy,numpoly}.isin(x, y)``

The definition of what those two features means needs defining before implementing.

``cumprod`` implementation

Same as #47, but quite a bit harder as even more bookeeping is needed. Might not even be possible without a element-by-element for-loop.

``prod`` implementation

Current implementation uses multiply repeatedly along axis. This can likely be
done much more efficiently with dedicated code.

Somewhat the same problem as issue #46, but a lot more book keeping.

Infinity loop in poly_divmod when high relative error is introduced

Issue in Chaospy: jonathf/chaospy#257

Problem with cross_truncate

I found a problem with cross_truncation. A concrete case is:
cross_truncate(numpy.array([[5,1,1,1,1]]), 10, 1.0)
We are expecting the answer to be true (the sum is 9 which is below 9), but the answer is False. I think this is due to numerical errors, and maybe the definition of the norm should be reformulated (to avoid divisions)?

Updating the PyPi package

PR #111 removed a bug with the newest version of numpy, but unfortunately has not yet found its way into the official repo.

Would you be willing to update the PyPi package?

Thanks!

Deprecation pkg_resources

This is a duplicate of jonathf/chaospy#408 for this package.
I can also propose a PR for this one if it is fine for you.

Function for gathering all lead coefficients

Get the coefficient in from of the polynomial term with the highest power.

Need to support ordering to deal with multivariate polynomial.