code: def reduce_axis(data, reducer, block_reducer, mapper=None, axis=None,
                  blen=None, storage=None, create='array', **kwargs):
      
  
      # setup
      storage = _util.get_storage(storage)
      blen = _util.get_blen_array(data, blen)
      length = len(data)
      # normalise axis arg
      if isinstance(axis, int):
          axis = (axis,)
  
      # deal with 'out' kwarg if supplied, can arise if a chunked array is
      # passed as an argument to numpy.sum(), see also
      # https://github.com/cggh/scikit-allel/issues/66
      kwarg_out = kwargs.pop('out', None)
      if kwarg_out is not None:
          raise ValueError('keyword argument "out" is not supported')
  
      if axis is None or 0 in axis:
          # two-step reduction
          out = None
          for i in range(0, length, blen):
              j = min(i+blen, length)
              block = data[i:j]
              if mapper:
                  block = mapper(block)
              res = reducer(block, axis=axis)
              if out is None:
                  out = res
              else:
                  out = block_reducer(out, res)
          if np.isscalar(out):
              return out
          elif len(out.shape) == 0:
              return out[()]
          else:
              return getattr(storage, create)(out, **kwargs)
  
      else:
          # first dimension is preserved, no need to reduce blocks
          out = None
          for i in range(0, length, blen):
              j = min(i+blen, length)
              block = data[i:j]
              if mapper:
                  block = mapper(block)
              r = reducer(block, axis=axis)
              if out is None:
                  out = getattr(storage, create)(r, expectedlen=length, **kwargs)
              else:
                  out.append(r)
          return out
docs: Apply an operation to `data` that reduces over one or more axes.