Conditional Plotting with Numpy and Pyplot

Question

I am trying to plot a function that is defined conditionally. Specifically: U(x) = (2**delta)/((D-d)**delta)*(D/2 - (x-x0))**delta , for abs(x-x0) less than D/2 and 0 otherwise.

But my problem is that I want to have x, x0 as numpy arrays, because this is how I use them in the rest of my actual code.

I have set-up the following example:

import numpy as np
import matplotlib.pyplot as plt
AD = 0.2
D = 0.4
delta = 8

def Parabolic(x, delta, D, AD):
    x0 = np.round(x)
    tempx = np.abs(x-x0)
    tempD = D/2*np.ones(len(x))
    if tempx<tempD:
        return ((2**delta)/(D-AD)**delta)*(D/2 - (x-x0))**delta
    else:
        return 0

figure = plt.figure(figsize=(10,8), dpi=72)  
xmin = -1.0
xmax = 1.0
X = np.linspace(xmin,xmax,1000)
plt.plot(X, Parabolic(X, delta=8, D=0.4, AD=0.2))

Obviously this example does not work, since the line tempx<tempD raises the error that a Truth-value of a list is ambiguous.

I searched around the documentation of numpy and found the function np.less(tempx, tempD). But if I replace tempx < tempD with np.less(tempx, tempD) it still doesn't work since once again I am asking for a truth value of an entire list. I understand that the problem is not with numpy, but with my inability to understand how to use the logical functions that numpy provides.

I am sorry if this answered some way in another post, I searched in this forum but could not find something else besides the curve() method. However I want to keep my numpy.array format for use in my actual codes. I would bet that the answer must be very simple, I just can't think of it.

On your error, see stackoverflow.com/a/30081652/901925 and its sidebar. — hpaulj
– hpaulj, Commented Apr 1, 2016 at 16:19

Xirtaminu · Accepted Answer · 2016-04-01 15:22:33Z

4

Try this which uses numpy logical arrays:

import numpy as np
import matplotlib.pyplot as plt
AD = 0.2
D = 0.4
delta = 8

def Parabolic(x, delta, D, AD):
    rtn_arr = np.zeros(len(x))
    x0 = np.round(x)
    tempx = np.abs(x-x0)
    tempD = D/2*np.ones(len(x))
    lgc_arr = tempx<tempD
    x_cut = x[lgc_arr]
    x0_cut = x0[lgc_arr]
    rtn_arr[lgc_arr] = ((2**delta)/(D-AD)**delta)*(D/2 - (x_cut-x0_cut))**delta
    return rtn_arr

figure = plt.figure(figsize=(10,8), dpi=72)
xmin = -1.0
xmax = 1.0
X = np.linspace(xmin,xmax,1000)
plt.plot(X, Parabolic(X, delta=8, D=0.4, AD=0.2))

answered Apr 1, 2016 at 15:22

Xirtaminu

4142 silver badges5 bronze badges

Sign up to request clarification or add additional context in comments.

3 Comments

George Datseris Over a year ago

Well, this is exactly what I was looking for. I managed to solve my problem using python lists and append, but of course it wasn't that good. Yours was what I was looking. I did not know that one could use x_cut = x[lgc_array] in order to get only the elements that correspond to true! One last thing: Could you explain me how the line lgc_arr = tempx<tempD worked? Does numpy arrays automatically change how the relation/comparison operators work and returns automatically a numpy array?

hpaulj Over a year ago

tempx<tempD is returning a array in your code as well. But you try to use that result in a Python if context. Using any numby boolean array in a context that expects a scalar boolean produces the ValueError.

Xirtaminu Over a year ago

@GeorgeDatseris Exactly as you say. The numpy arrays redefine the conditional operators (which can be done for any custom class docs.python.org/2/library/stdtypes.html) so that using >, < etc. on numpy arrays returns numpy arrays. Then tempx<tempD is an array of boolean values, which can be used to do array slicing with x[lgc_array].

B. M. · Accepted Answer · 2016-04-01 21:11:59Z

3

Parabolic must be a ufunc, so you can't put python test in your code.

a simple workaround is :

def Parabolic(x, delta, D, AD):
    x0 = np.round(x)
    tempx = np.abs(x-x0)
    tempD = D/2*np.ones(len(x))
    u=(((2**delta)/(D-AD)**delta)*(D/2 - (x-x0))**delta)
    u[tempx>=tempD]=0
    return u

or to avoid unnecessary computations:

def Parabolic2(x, delta, D, AD):
    x0 = np.round(x)
    tempx = np.abs(x-x0)
    tempD = D/2*np.ones(len(x))
    u= zeros_like(x)
    valid=tempx<tempD
    u[valid]=(((2**delta)/(D-AD)**delta)*(D/2 - (x-x0)[valid])**delta)
    return u

The second is slighty faster :

In [141]: %timeit Parabolic(x,8,.4,.2)
1000 loops, best of 3: 310 µs per loop

In [142]: %timeit Parabolic2(x,8,.4,.2)
1000 loops, best of 3: 218 µs per loop

edited Apr 1, 2016 at 21:11

answered Apr 1, 2016 at 15:31

B. M.

18.7k2 gold badges40 silver badges56 bronze badges

3 Comments

George Datseris Over a year ago

This is also a nice workaround, and quite simple code-wise. However it does mean that I have to calculate the function for all values, small and big ones. And this will take a lot of time in the normal code. I intend to use the first answer for both very high and very small values of the u function, mainly to avoid computations.

Xirtaminu Over a year ago

I don't think the second version will work as the length of u[tempx<tempD] will be less than the array returned on the right side of this equation.

George Datseris Over a year ago

I tried to run the code and it is true that it doesn't complete for the reason you said.

Collectives™ on Stack Overflow

Conditional Plotting with Numpy and Pyplot

2 Answers 2

3 Comments

3 Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

3 Comments

3 Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related