code: def per_base(x, windows, is_accessible=None, fill=np.nan):
      
  
      # calculate window sizes
      if is_accessible is None:
          # N.B., window stops are included
          n_bases = np.diff(windows, axis=1).reshape(-1) + 1
      else:
          n_bases = np.array([np.count_nonzero(is_accessible[i-1:j])
                              for i, j in windows])
  
      # deal with multidimensional x
      if x.ndim == 1:
          pass
      elif x.ndim == 2:
          n_bases = n_bases[:, None]
      else:
          raise NotImplementedError('only arrays of 1 or 2 dimensions supported')
  
      # calculate density per-base
      with ignore_invalid():
          y = np.where(n_bases > 0, x / n_bases, fill)
  
      # restore to 1-dimensional
      if n_bases.ndim > 1:
          n_bases = n_bases.reshape(-1)
  
      return y, n_bases
docs: Calculate the per-base value of a windowed statistic.
  
      Parameters
      ----------
  
      x : array_like, shape (n_windows,)
          The statistic to average per-base.
      windows : array_like, int, shape (n_windows, 2)
          The windows used, as an array of (window_start, window_stop)
          positions using 1-based coordinates.
      is_accessible : array_like, bool, shape (len(contig),), optional
          Boolean array indicating accessibility status for all positions in the
          chromosome/contig.
      fill : object, optional
          Use this value where there are no accessible bases in a window.
  
      Returns
      -------
  
      y : ndarray, float, shape (n_windows,)
          The input array divided by the number of (accessible) bases in each
          window.
      n_bases : ndarray, int, shape (n_windows,)
          The number of (accessible) bases in each window