docs: 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