2

I'm setting a numpy array with a power-law equation. The problem is that part of my domain tries to do numpy.power(x, n) when x is negative and n is not an integer. In this part of the domain I want the value to be 0.0. Below is a code that has the correct behavior, but is there a more Pythonic way to do this?

# note mesh.x is a numpy array of length nx
myValues = npy.zeros((nx))
para = [5.8780046, 0.714285714, 2.819250868] 
for j in range(nx):
    if mesh.x[j] > para[1]:
        myValues[j] = para[0]*npy.power(mesh.x[j]-para[1],para[2])

    else:
        myValues[j] = 0.0
Improve this question

3 Answers 3

Reset to default
1

Is "numpythonic" a word? It should be a word. The following is really neither pythonic nor unpythonic, but it is much more efficient than using a for loop, and close(r) to the way Travis would probably do it:

import numpy
mesh_x = numpy.array([0.5,1.0,1.5])
myValues = numpy.zeros_like( mesh_x )
para = [5.8780046, 0.714285714, 2.819250868] 
mask = mesh_x > para[1]
myValues[mask] = para[0] * numpy.power(mesh_x[mask] - para[1], para[2])
print(myValues)

For very large problems you would probably want to avoid creating temporary arrays:

mask = mesh.x > para[1]
myValues[mask] = mesh.x[mask]
myValues[mask] -= para[1]
myValues[mask] **= para[2]
myValues[mask] *= para[0]
Improve this answer
Sign up to request clarification or add additional context in comments.

3 Comments

The first clips throws an exception. Traceback (most recent call last): File "C:\Users\choutman\Documents\confocal beads\Modelling\Carl_1D_coupled_time.py", line 85, in <module> myValues[mask] = para[0] * npy.power(mesh.x[mask] - para[1], para[2]) {{clip}}File "C:\Python27\lib\site-packages\fipy\variables\unaryOperatorVariable.py", line 67, in calcValue return self.op(self.var[0].value) File "C:\Python27\lib\site-packages\fipy\variables\variable.py", line 1572, in <lambda> return self._UnaryOperatorVariable(lambda a: a[index], IndexError: unsupported iterator index
I cannot replicate your problem. I have edited my answer to include a complete example that runs for me. You might get problems if your mesh.x is not a standard numpy.array, but rather a matrix or other derived class - could this be the case? If so, maybe say mask=(numpy.asarray(mesh.x)>para[1]) to be safe.
Yes, exactly. mesh.x is a class derived from FiPy. When I cast it as a numpy array as you suggest, it works fine. Thanks.
1

Here's one approach with np.where to choose values between the power calculations and 0 -

import numpy as np
np.where(mesh.x>para[1],para[0]*np.power(mesh.x-para[1],para[2]),0)

Explanation :

  1. np.where(mask,A,B) chooses elements from A or B depending on mask elements. So, in our case it is mesh.x>para[1] when doing a vectorized comparison for all mesh.x elements in one go.
  2. para[0]*np.power(mesh.x-para[1],para[2]) gives us the elements that are to be chosen in case a mask element is True. Else, we choose 0, which is the third argument to np.where.
Improve this answer

7 Comments

+1 for teaching me about this usage of numpy.where that I wasn't aware of, but back to 0 for the fact that the subtraction, exponentiation and multiplication happen for all elements of mesh.x which are then selected post-hoc. Would this not make the result a complex array, because of the cases that violate the condition?
uhh, apparently not, if mesh.x isn't complex to begin with: instead it returns nan in the bad locations. Again I learn something. But you do get RuntimeWarning: invalid value encountered in power
@jez I think those bad locations would be covered up by the else part of 0 filling, right?
yes, they are covered up but the calculation still happens (wasting CPU time and issuing the warning) for the bad values. That's what I mean by "selected post-hoc".
I first tried the np.where and got this: C:\Users\choutman\Documents\confocal beads\Modelling\Carl_1D_coupled_time.py:84: RuntimeWarning: invalid value encountered in power npy.where(mesh.x>para[1],para[0]*npy.power(mesh.x-para[1],para[2]),0) C:\Python27\lib\site-packages\fipy\variables\binaryOperatorVariable.py:81: RuntimeWarning: invalid value encountered in power return self.op(self.var[0].value, val1)
|
0

More of an explanation of the answers given by @jez and @Divakar with simple examples than an answer itself. They both rely on some form of boolean indexing.

>>> 
>>> a
array([[-4.5, -3.5, -2.5],
       [-1.5, -0.5,  0.5],
       [ 1.5,  2.5,  3.5]])
>>> n = 2.2
>>> a ** n
array([[         nan,          nan,          nan],
       [         nan,          nan,   0.21763764],
       [  2.44006149,   7.50702771,  15.73800567]])

np.where is made for this it selects one of two values based on a boolean array.

>>> np.where(np.isnan(a**n), 0, a**n)
array([[  0.        ,   0.        ,   0.        ],
       [  0.        ,   0.        ,   0.21763764],
       [  2.44006149,   7.50702771,  15.73800567]])
>>> 
>>> b = np.where(a < 0, 0, a)
>>> b
array([[ 0. ,  0. ,  0. ],
       [ 0. ,  0. ,  0.5],
       [ 1.5,  2.5,  3.5]])
>>> b **n
array([[  0.        ,   0.        ,   0.        ],
       [  0.        ,   0.        ,   0.21763764],
       [  2.44006149,   7.50702771,  15.73800567]])

Use of boolean indexing on the left-hand-side and the right-hand-side. This is similar to np.where

>>> 
>>> a[a >= 0] = a[a >= 0] ** n
>>> a
array([[ -4.5       ,  -3.5       ,  -2.5       ],
       [ -1.5       ,  -0.5       ,   0.21763764],
       [  2.44006149,   7.50702771,  15.73800567]])
>>> a[a < 0] = 0
>>> a
array([[  0.        ,   0.        ,   0.        ],
       [  0.        ,   0.        ,   0.21763764],
       [  2.44006149,   7.50702771,  15.73800567]])
>>>
Improve this answer

Comments

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.