Whoops, I forgot that I asked about this on the mailing list two years ago: https://mail.scipy.org/pipermail/numpy-discussion/2015-February/072183.html

(nabble: http://numpy-discussion.10968.n7.nabble.com/Any-interest-in-a-heaviside-ufunc-td39641.html)

Implementation seems sound to me, but I share some of the sentiments in the mailing list about whether numpy is the right place for this

I think this would be a reasonable addition. The other versions of the Heaviside function (unitstep functions) are easy to implement, but this version with .5 at zero is not so straight forward. Of the various signal waveforms, rectangular, hat, etc., this is probably the only one that belongs in numpy.

How hard would it be to allow default values in ufuncs? A default of 0.5 seems useful to me?

A default of 0.5 seems useful to me?

Perhaps, but I think this is likely to catch people out, and so being explicit would be better.

Defaults on ufuncs (and also, gufuncs) is something that might be useful elsewhere though

Also, perhaps worth linking to the wikipedia page about the choice of H(0) in the documentation?

I note that the .5 value also comes from Fourier series of piecewise smooth functions with discrete discontinuities, which you can see, along with the Gibbs phenomenon, from convolution with the Dirichlet kernels,

@seberg I think adding arguments to ufuncs will be a problem.

I added a "References" section to the docstring, with a link to the Wikipedia article.

Adding default values for arguments doesn't appear to be possible with the existing ufunc-generating machinery. (A big +1 for the enhancement request, though.)

How hard would it be to allow default values in ufuncs? A default of 0.5 seems useful to me?

I think dealing with multiple types would be a problem in general. With scalar defaults the same kind casting rules would probably work in this case where the only valid types are floats. Although it would still be a problem for scalar values of x.

OK, let's put this in. Thanks @WarrenWeckesser .

@WarrenWeckesser There is probably no reason not to allow integer types here, that way if h0 is set to 0 or 1 there is a nice clipping function.

@charris, do you mean something like heaviside([-1, 0, 1], 1) should return the integer array array([ 0, 1, 1])?

This needs to be listed in the documentation somewhere too, presumably routines.umath.rst - otherwise it won't be discoverable

@eric-wieser Looks like it should be in docs/source/reference/ufuncs.rst and docs/source/reference/routines.math.rst.

I'll get started on a new pull request to update the docs and handle integers.

@charris (or anyone else): Can you recommend an existing ufunc implementation similar to heaviside that I can study/copy so that it handles integers and floats? I'm running into problems getting the TD() declarations correct in generate_umath.py. I have two (or three) C functions, one for floating point values and one for integers (or two for integers: one each for signed and unsigned integers).

I'll keep poking at this, but some pointers might save me some time.

Might my changes in #8774 help? If you name your functions the right thing, I don't think you need any real work in generate_umath.py

Thanks. I don't plan on implementing the code for handling generic Python objects. I've been using the following in npy_math_internal.h.src:

/**begin repeat
 * #type = npy_float, npy_double, npy_longdouble#
 * #c = f,d,g#
 */

@type@ npy_heaviside@c@(@type@ x, @type@ h0)
{
    if (npy_isnan(x)) {
        return (@type@) NPY_NAN;
    }
    else if (x == 0) {
        return h0;
    }
    else if (x < 0) {
        return (@type@) 0.0;
    }
    else {
        return (@type@) 1.0;
    }
}

/**end repeat**/


/**begin repeat
 * #type = npy_byte, npy_short, npy_int, npy_long, npy_longlong#
 * #c = b,h,i,l,q#
 */

@type@ npy_heaviside@c@(@type@ x, @type@ h0)
{
    if (x == 0) {
        return h0;
    }
    else if (x < 0) {
        return (@type@) 0;
    }
    else {
        return (@type@) 1;
    }
}

/**end repeat**/


/**begin repeat
 * #type = npy_ubyte, npy_ushort, npy_uint, npy_ulong, npy_ulonglong#
 * #c = B,H,I,L,Q#
 */

@type@ npy_heaviside@c@(@type@ x, @type@ h0)
{
    if (x == 0) {
        return h0;
    }
    else {
        return (@type@) 1;
    }
}

/**end repeat**/

(I have no problem dropping the specialized code for the unsigned integers if that seems like overkill.)

What is the appropriate way to add this to defdict in generate_umath.py? I haven't figured out the appropriate instantiation of a TD() (or TypeDescription()s) to add to the Ufunc() object.

My point was that exposing a function called @TYPE@_heaviside for each integer type as I do here means that a simple TD(ints) will suffice

Good timing! After looking at your pull request some more, I was just coming to that realization. That code looks very boilerplate-ish. I guess I expected code at that level to be generated by the python script, especially after implementing the "float only" code without making such changes.

Yeah, to be honest I wasn't aware that the way you're doing it was possible either, so I guess we're even there. But the way I'm doing it is consistent with how np.maximum does it, so I'm not reinventing a wheel - just perhaps using one that was supposed to have been removed.

@charris wrote

There is probably no reason not to allow integer types here, that way if h0 is set to 0 or 1 there is a nice clipping function.

How important do folks think it is that the function return an integer array if the inputs are all integers? Currently it always returns floats. I started looking into this, but I hit a steeper section of the ufunc implementation learning curve, and then got caught up in other things.

It could be changed after the next release, but that means there would be a change in behavior of, for example, heaviside([-1, 0, 1], 1) when the code started handling integers differently. What is the cutoff date for 0.13 code?

What is the cutoff date for 0.13 code?

Four weeks ago ;) Seriously, all I want is to finish up __array_ufunc__ and maybe merge some left over small bits. Oh, and figure out what is going on with the scipy errors.

On my machine, compiling master, heaviside is an order of magnitude slower than a naive cast of x > 0:

In [1]: x = np.random.randn(100000)

In [2]: %timeit np.heaviside(x, 0)
1.31 ms ± 58.3 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

In [3]: %timeit np.where(x > 0, 1.0, 0.0)
658 µs ± 74.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

In [4]: %timeit (x>0).astype(x.dtype)
172 µs ± 34.8 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

Right now, it would seemingly be better to have defined heaviside as a normal function.

Possible improvements:

• Write a direct loop, rather than deferring to function pointers
• add a fast-path for the almost certain case that x0 is a scalar
• Further special case the values of 0, 0.5, and 1.0 in those cases

One more possible improvement to add to the list:

• Make the return type agree with the input types (i.e. don't always return floats).

I wonder if this shouldn't be in scipy.signals

Moreover it could make sense to have also a parameter for "delayed" heaviside function.

PS: sorry for noise... about delayed heaviside function...

t = np.arange(-1.0, 10.0, 0.01)
u = np.heaviside(t - 4, 1.0)
plt.plot(t, u)

does the trick (for a 4s delay)