Say we want a function like lerp() to work for t a scalar or range, and p0 p1 scalars or arrays of any shape:

import numpy as np

def lerp( t, p0, p1 ):
    return (1 - t) * p0 + t * p1

p0 = np.array([ 1, 2, 3 ])
p1 = np.array([ 4, 5, 6 ])
T = np.linspace( .1, .5, 5 )

    # We want lerp( T, p0, p1 ) to be [lerp(.1), lerp(.2) ... lerp(.5)], but:
print lerp( T, p0, p1 )
    # ValueError: shape mismatch: objects cannot be broadcast to a single shape

    # Worse,
print lerp( np.linspace( .1, .5, 3 ), p0, p1 )
    # => a nonsense result, with no warning

"Broadcasting" is the part of numpy that expands e.g. 2*p0 to array([ 2, 4, 6 ]), and p0*p1 to
array([ 4, 10, 18 ]). However, broadcasting leads to naive code like the above working for some arguments but not for others. A way to restrict broadcasting in such cases, from Josef, is to make t a column vector:

def lerp( t, p0, p1 ):
    if not np.isscalar(t) and not np.isscalar(p0): 
        t = t[:,None]  # make t a column vec, to restrict broadcasting
    return (1 - t) * p0 + t * p1

print lerp( T, p0, p1 )  # [lerp(.1), lerp(.2) ... lerp(.5)]
print lerp( np.linspace( .1, .5, 3 ), p0, p1 )  # [lerp(.1), lerp(.3), lerp(.5)]

A bug, though: if p0 p1 are column vectors of the same size as t, t*p will broadcast to the elementwise product, which is not what you want.

Links:
EricsBroadcastingDoc
.../site-packages/numpy/doc/broadcasting.py
Haeberli and Voorhies, Image Processing By Interpolation and Extrapolation
-- nothing to do with broadcasting, just a lovely paper on Lerp