[ENH]: Expand shape compatibility of x and y in plot

[jsdodge]

Problem

When x is 1D and y is 2D, the expression plot(x, y) behaves differently in Matplotlib and MATLAB. Matplotlib requires that y.shape[0] == x.size and always plots the columns of y against x. MATLAB will plot by y by columns if its first dimension is compatible with x, but if it is not, and the second dimension of y is compatible with x, then it will plot y by rows.

The following code works in MATLAB, for example:

x = linspace(0, 2*pi);    % Size: [1 x 100]
k = (0:2)';               % Size: [3 x 1]
y = sin(k*x);             % Size: [3 x 100]
plot(x, y)                % Second dimensions are the same
plot(x, y')               % Second dimension of x is the same as first dimension of y'

In the similar code below, Matplotlib returns an error for plot(x, y), but not for plot(x, y.T):

x = np.linspace(0, 2 * np.pi)    # Shape: (50,)
k = np.arange(3).                # Shape: (3,)
y = np.sin(np.outer(k, x))       # Shape: (3, 50)
plt.plot(x, y)                   # Returns ValueError
plt.plot(x, y.T)                 # y.T has shape (50, 3), compatible with x

The NumPy broadcasting rules and FFT routines are organized around row operations, which forces a choice between rows and columns—see here for an example.

Proposed solution

When x is 1D and y is 2D, plot y by columns when y.shape[0] == x.size, but also try to plot y by rows when y.shape[0] != x.size.

[timhoffm]

Thanks for the propsal. However, I believe supporting this would be an API flaw. "y columns represent the datasets" is a clear definition. Adding "but if this does not match the number of points in x, we try to use the rows as datasets" results in magical behavior. Note also that you can plot only y in which case plot(y) might result in different datasets that plot(x, y) 😨.

OTOH, you should know your data and can easily make the row into columns by explicitly transposing plot(x, y.T) - "explicit is better than implicit".

[jsdodge]

Thanks for your reply, @timhoffm. I see your point, but assign more value to the convenience of being able to ignore the row/column distinction when I just want to make a simple plot of a few lines, and consider support for that to be "magical" in a positive way. The use case of plotting a 2D y versus a 1D x is so common that I think it would be worth making an exception to the "explicit is better than implicit" guideline in favor of the "practicality beats purity" guideline. The MATLAB developers have implemented this behavior for a reason.

Regarding the potential inconsistency between plot(x, y) and plot(y), I don't see how this would be a problem in practice. If y has shape (r, c), then plot(y) will produce c lines of r data points. If that is what the user wants, then they would also use an x that has shape (r,) with plot(x, y). The only ambiguity that can appear is when y is a square 2D array that the user wants to plot by rows with plot(y.T), and in that case the conflict would be between plot(x, y) and plot(y.T). In my experience, most of the time the number of lines is much smaller than the number of points in a line, so even this case seems unlikely to occur with any frequency.

For signal processing applications, the current behavior also requires users to transpose processed 2D data almost by default, since the NumPy FFT routines are applied along the last index of an ND array by default and the NumPy broadcasting rules also privilege the last index. The code below, taken from my Discourse post, provides an example. MATLAB applies the FFT along the first index, which is better integrated with the column orientation of the MATLAB-inspired Matplotlib plot function.

ns = 256
nw = 10
dt = 10 / ns
t = np.arange(-ns * dt / 2, ns * dt / 2, dt)
x = gausspulse(t, fc=1)

rng = default_rng(0)

# Broadcasting matches rightmost indices
xm = x + 0.1 * np.abs(x) * rng.standard_normal((nw, ns))

# Default axis for FFT is -1
f = rfftfreq(ns, t[1] - t[0])
xm_f = rfft(xm)

# Need to take the transpose of xm and xm_f to convert rows to columns
_, axs = plt.subplots(2, 1)
axs[0].plot(t, xm.T)
axs[1].plot(f, np.abs(xm_f.T) ** 2)
plt.show()

[anntzer]

Rght now, numpy only lets you broadcast against specific dimensions, i.e.

• if x has shape (10,) and y has shape (20, 10), then the expression x+y is valid and has shape (20, 10), but
• if x has shape (20,) and y has shape (20, 10), then the expression x+y is not valid, even though you could argue that it should have shape (20, 10).

If you can convince the numpy devs to change the rule above and introduce "size-matching" to allow size(20,) + size(20, 10), then we could reconsider the behavior on our side. If numpy doesn't change that rule (which seems exceedingly unlikely, to be honest), this will realistically not happen.