def smoothness(y, r=20, stepsize=1, domain=None, unsigned=False):
    """
    This fucntion calculates a smoothness measure on the input 1D data.
    The concept is similar to that decribed here: http://www.goldensoftware.com/variogramTutorial.pdf

    USAGE
    ---
    u,x = smoothness(t,y,r=4)

    INPUTS
    ---
    t,          Domain points of data set
    y,          Range of data set
    r=4,        Radius of which to consider variations (derivates)
    stepsize=1, The average will be considered every stepsize points. Increasing
                this from its default of 1 can be useful when processing large
                data sets; however, when stepsize is not 1, the length of the
                output will differ from that of the inputs.

    OUTPUTS
    ---
    u,       Sub-domain which is one-to-one with smoothness measure
    x,       Smoothness measure -- the data, y, is smooth when |x| is approx. 1
             NOTE that x=-1 is smooth and decreasing while x=1 is smooth and increasing

    """
    from numpy import arange, var, std, polyfit, poly1d, mean, diff, zeros_like, array
    from scipy.interpolate import InterpolatedUnivariateSpline as spline
    if domain is None:
        domain = list(range(0, len(y)))
    x, u = [], []
    for k in arange(0, len(y), stepsize):
        a = max(0, k - r)
        b = min(len(y), k + r) - 1
        D = (y[b] - y[a]) / (b - a)
        if unsigned:
            D = abs(D)
        d = abs(mean(diff(y[a:b])))
        x.append(D / d if d != 0 else 0)
        u.append((domain[a] + domain[b]) / 2)
    x = array(x)
    if stepsize > 1:
        x = spline(u, x, k=1)(domain)
    return x
