docs: def get_gradient_y(shape, py):
      """Calculate the gradient in the y direction to the line at py
  
      The y gradient operator is a block matrix, where each block is the size of the image width.
      The matrix itself is made up of (img_height x img_height) blocks, most of which are all zeros.
      """
      import scipy.sparse
  
      height, width = shape
      rows = []
      empty = scipy.sparse.dia_matrix((width, width))
      identity = scipy.sparse.identity(width)
  
      # Create the blocks by row, beginning with blocks leading up to the peak row from the top
      for n in range(py):
          row = [empty]*n
          row += [-identity, identity]
          row += [empty]*(height-n-2)
          rows.append(row)
      # Set all elements in the peak row to zero
      rows.append([empty]*height)
      # Create the blocks for the rows leading up to the peak row from the bottom
      for n in range(height-py-1):
          row = [empty]*(py+n)
          row += [identity, -identity]
          row += [empty]*(height-py-n-2)
          rows.append(row)
      return scipy.sparse.bmat(rows)