def regularized_function(x,y,func,bins=None,range=None):
    """Compute func() over data aggregated in bins.

    (x,y) --> (x', func(Y'))  with Y' = {y: y(x) where x in x' bin}

    First the data is collected in bins x' along x and then func is applied to
    all data points Y' that have been collected in the bin.

    :Arguments:
       x
          abscissa values (for binning)
       y
          ordinate values (func is applied)
       func
          a numpy ufunc that takes one argument, func(Y')
       bins
          number or array
       range
          limits (used with number of bins)

    :Returns:
       F,edges
          function and edges (midpoints = 0.5*(edges[:-1]+edges[1:]))
    """
    _x = numpy.asarray(x)
    _y = numpy.asarray(y)

    # setup of bins from numpy.histogram
    if (range is not None):
        mn, mx = range
        if (mn > mx):
            raise ValueError('max must be larger than min in range parameter.')

    if not numpy.iterable(bins):
        if range is None:
            range = (_x.min(), _x.max())
        mn, mx = [float(mi) for mi in range]
        if mn == mx:
            mn -= 0.5
            mx += 0.5
        bins = numpy.linspace(mn, mx, bins+1, endpoint=True)
    else:
        bins = numpy.asarray(bins)
        if (numpy.diff(bins) < 0).any():
            raise ValueError('bins must increase monotonically.')

    sorting_index = numpy.argsort(_x)
    sx = _x[sorting_index]
    sy = _y[sorting_index]

    # boundaries in SORTED data that demarcate bins; position in bin_index is the bin number
    bin_index = numpy.r_[sx.searchsorted(bins[:-1], 'left'),
                         sx.searchsorted(bins[-1], 'right')]

    # naive implementation: apply operator to each chunk = sy[start:stop] separately
    #
    # It's not clear to me how one could effectively block this procedure (cf
    # block = 65536 in numpy.histogram) because there does not seem to be a
    # general way to combine the chunks for different blocks, just think of
    # func=median
    F = numpy.zeros(len(bins)-1)  # final function
    F[:] = [func(sy[start:stop]) for start,stop in izip(bin_index[:-1],bin_index[1:])]
    return F,bins
Compute func() over data aggregated in bins.

    (x,y) --> (x', func(Y'))  with Y' = {y: y(x) where x in x' bin}

    First the data is collected in bins x' along x and then func is applied to
    all data points Y' that have been collected in the bin.

    :Arguments:
       x
          abscissa values (for binning)
       y
          ordinate values (func is applied)
       func
          a numpy ufunc that takes one argument, func(Y')
       bins
          number or array
       range
          limits (used with number of bins)

    :Returns:
       F,edges
          function and edges (midpoints = 0.5*(edges[:-1]+edges[1:]))