def piecewise_integrate(x, y, a, b):
    """
    Integrate a piecewise constant function.

    Parameters
    ----------
    x : 1D numpy array
        Sorted distinct values. The function we integrate changes at these x
        values.

    y : 1D numpy array
        The function we are integrating specified for each value in x.
        Specifically, the value from [x[0], x[1]) is y[0], etc.

    a : float
        This is where we start integrating from. The code here actually
        requires a to be equal to x[0].

    b : float
        This is where we stop integrating.

    Returns
    -------
    output : float
        The computed integral.
    """
    assert x[0] == a
    assert x[-1] <= b
    output = 0.
    num_x = len(x)
    if x[-1] == b:
        for idx in range(num_x - 1):
            output += y[idx] * (x[idx+1] - x[idx])
    else:
        for idx in range(num_x):
            if idx < num_x - 1:
                output += y[idx] * (x[idx+1] - x[idx])
            else:
                output += y[idx] * (b - x[idx])
    return output