How to raise arrays with negative values to fractional power in Python?

Question

I have an array with negative values that has to be raised to fractional power in Python. I need to obtain the real part of the complex number array generated by the operation.

MWE

from __future__ import division
import numpy as np
a = -10
b = 2.5
n = 0.88
x = np.arange(5, 11, 1)
y = (a / (x - b)) ** (1 / n)

I am using Python v2.7.6.

Answer 1

The issue is that NumPy does not promote float or integer dtypes to complex dtypes for this calculation.

You have a float array base and a float exponent, so NumPy tries to compute the results using the "put two float dtype objects in, get a float dtype object out" loop. Negative values trigger a warning and returns an array of null values.

** end up using the same code as np.power when one of the operands is an array. You can see all of the low-level loops that can be used below. Note that you always get back an object with the same dtype as your input dtypes:

>>> np.power.types
['bb->b', # char
 'BB->B', # unsigned char
 'hh->h', # short
  ...     
 'dd->d', # compatible Python float
 'gg->g', # compatible: C long float
 'FF->F',
 'DD->D', # compatible: Python complex
 'GG->G',
 'OO->O']

We want the calculation to run with the 'DD->D' loop!

The solution, as pointed out by others on this page, is to make sure that either the base or the exponent has a complex dtype. This forces NumPy to promote any lesser numeric dtypes to a complex dtype and the computation uses the "put two complex dtype objects in, get a complex dtype object out" loop:

>>> a = -10 + 0j # ensures that a/(x - b) will have complex dtype
>>> ((a / (x - b)) ** (1 / n))
array([-4.39566725-2.00743397j, -2.99895689-1.36957772j,
       -2.25394034-1.02934006j, -1.79435400-0.81945401j,
       -1.48410349-0.67776735j, -1.26136729-0.57604714j])

If you just want the real parts use ((a / (x - b)) ** (1 / n)).real.

Answer 2

As you are using numpy, you could use np.power as well like this:

from __future__ import division
import numpy as np
a = -10+0j
b = 2.5
n = 0.88
x = np.arange(5, 11, 1)
y = np.power((a / (x - b)),(1 / n))

However, actually ** of an numpy array is just syntactic sugar for ndarray.pow() and will result in the same code being executed as np.power. See @acjr's comment.