docs: def apply_once(func, arr, axes, keepdims=True):
      """
      Similar to `numpy.apply_over_axes`, except this performs the operation over
      a flattened version of all the axes, meaning that the function will only be
      called once. This only makes a difference for non-linear functions.
  
      Parameters
      ----------
      func : callback
          Function that operates well on Numpy arrays and returns a single value
          of compatible dtype.
      arr : ndarray
          Array to do operation over.
      axes : int or iterable
          Specifies the axes to perform the operation. Only one call will be made
          to `func`, with all values flattened.
      keepdims : bool
          By default, this is True, so the collapsed dimensions remain with
          length 1. This is simlar to `numpy.apply_over_axes` in that regard.  If
          this is set to False, the dimensions are removed, just like when using
          for instance `numpy.sum` over a single axis. Note that this is safer
          than subsequently calling squeeze, since this option will preserve
          length-1 dimensions that were not operated on.
  
      Examples
      --------
      >>> import deepdish as dd
      >>> import numpy as np
      >>> rs = np.random.RandomState(0)
      >>> x = rs.uniform(size=(10, 3, 3))
  
      Image that you have ten 3x3 images and you want to calculate each image's
      intensity standard deviation:
  
      >>> np.apply_over_axes(np.std, x, [1, 2]).ravel()
      array([ 0.06056838,  0.08230712,  0.08135083,  0.09938963,  0.08533604,
              0.07830725,  0.066148  ,  0.07983019,  0.08134123,  0.01839635])
  
      This is the same as ``x.std(1).std(1)``, which is not the standard
      deviation of all 9 pixels together. To fix this we can flatten the pixels
      and try again:
  
      >>> x.reshape(10, 9).std(axis=1)
      array([ 0.17648981,  0.32849108,  0.29409526,  0.25547501,  0.23649064,
              0.26928468,  0.20081239,  0.33052397,  0.29950855,  0.26535717])
  
      This is exactly what this function does for you:
  
      >>> dd.apply_once(np.std, x, [1, 2], keepdims=False)
      array([ 0.17648981,  0.32849108,  0.29409526,  0.25547501,  0.23649064,
              0.26928468,  0.20081239,  0.33052397,  0.29950855,  0.26535717])
      """
  
      all_axes = np.arange(arr.ndim)
      if isinstance(axes, int):
          axes = {axes}
      else:
          axes = set(axis % arr.ndim for axis in axes)
  
      principal_axis = min(axes)
  
      for i, axis in enumerate(axes):
          axis0 = principal_axis + i
          if axis != axis0:
              all_axes[axis0], all_axes[axis] = all_axes[axis], all_axes[axis0]
  
      transposed_arr = arr.transpose(all_axes)
  
      new_shape = []
      new_shape_keepdims = []
      for axis, dim in enumerate(arr.shape):
          if axis == principal_axis:
              new_shape.append(-1)
          elif axis not in axes:
              new_shape.append(dim)
  
          if axis in axes:
              new_shape_keepdims.append(1)
          else:
              new_shape_keepdims.append(dim)
  
      collapsed = np.apply_along_axis(func,
                                      principal_axis,
                                      transposed_arr.reshape(new_shape))
  
      if keepdims:
          return collapsed.reshape(new_shape_keepdims)
      else:
          return collapsed