CODEWITHSUNDEEP
pythonarraysnumpypolynomials
Sebastian Mendez
Sebastian Mendez

Reputation: 2981

Numpy Problems with Arrays of poly1d Objects

I'd like to first start this out with the fact that it is possible, in numpy, to create an array of poly1d objects:

random_poly = np.frompyfunc(lambda i, j: np.poly1d(np.random.randint(1, 4, 3)), 2, 1)
def random_poly_array(shape):    
    return np.fromfunction(random_poly, shape)

a1 = random_poly_array((3,3))

This works just fine, and we can even multiply matrices made from this form using np.dot:

a2 = random_poly_array((3,3))
a1_x_a2 = np.dot(a1, a2)

However, most other methods fail to work. For example, you can't take a list of certain poly1d objects and convert it into an array:

np.array([np.poly1d([1,2,3]), np.poly1d([1,2,3])])

As that will raise ValueError: cannot copy sequence with size 2 to array axis with dimension 3. To add to the confusion,

np.array([np.poly1d([1,2]), np.poly1d([1,2])])

will not raise an error, but instead create a 2x2 array of just 2's. Adding dtype=object has no affect, and numpy will still try to convert the poly1d objects to arrays.

The reason why this is problematic is that one cannot take an array of dimension d and convert it to an array of poly1d objects of dimension d-1. I would have expected

arr = np.arange(1, 10).reshape(3,3)
np.apply_along_axis(np.poly1d, 0, arr)

To return an array of poly1d objects, but instead it returns an unalterated array. Even worse, if arr=np.arange(9).reshape(3,3), it will throw an error, as the first poly1d object created will have a length of 2 instead of 3 due to the zero coefficient. Thus, my question is this: is there a feasible method to create poly1d arrays in numpy? If not, why not?

Upvotes: 1

Views: 908

Answers (1)

Sebastian Mendez
Sebastian Mendez

Reputation: 2981

Using the concept of None forcing numpy to not broadcast an object into an array, something brought to my attention by Paul Panzer, I created a function which will transform the last axis into a poly1d object:

def array_to_poly(arr):
    return np.apply_along_axis(lambda poly: [None, np.poly1d(poly)], -1, arr)[..., 1]

However, if we're okay with abusing more than one system in a single function, we can make it apply over arbitrary axes:

def array_to_poly(arr, axis=-1):
    temp_arr = np.apply_along_axis(lambda poly: [None, np.poly1d(poly)], axis, arr)
    n = temp_arr.ndim
    s = [slice(None) if i != axis%n else 1 for i in range(n)]
    return temp_arr[s]

Testing it with arr = np.arange(1, 25).reshape(2,3,4), we obtain:

In [ ]: array_to_poly(arr, 0)
Out[ ]: 
array([[poly1d([ 1, 13]), poly1d([ 2, 14]), poly1d([ 3, 15]),
        poly1d([ 4, 16])],
       [poly1d([ 5, 17]), poly1d([ 6, 18]), poly1d([ 7, 19]),
        poly1d([ 8, 20])],
       [poly1d([ 9, 21]), poly1d([10, 22]), poly1d([11, 23]),
        poly1d([12, 24])]], dtype=object)

In [ ]: array_to_poly(arr, 1)
Out[ ]: 
array([[poly1d([1, 5, 9]), poly1d([ 2,  6, 10]), poly1d([ 3,  7, 11]),
        poly1d([ 4,  8, 12])],
       [poly1d([13, 17, 21]), poly1d([14, 18, 22]), poly1d([15, 19, 23]),
        poly1d([16, 20, 24])]], dtype=object)

In [ ]: array_to_poly(arr, 2)
Out[ ]: 
array([[poly1d([1, 2, 3, 4]), poly1d([5, 6, 7, 8]),
        poly1d([ 9, 10, 11, 12])],
       [poly1d([13, 14, 15, 16]), poly1d([17, 18, 19, 20]),
        poly1d([21, 22, 23, 24])]], dtype=object)

as expected.

Upvotes: 1

Related Questions