docs: def bound_perturbed_gmres(pseudo, p, epsilon, deltas):
      '''Compute GMRES perturbation bound based on pseudospectrum
  
      Computes the GMRES bound from [SifEM13]_.
      '''
      if not numpy.all(numpy.array(deltas) > epsilon):
          raise ArgumentError('all deltas have to be greater than epsilon')
  
      bound = []
      for delta in deltas:
          # get boundary paths
          paths = pseudo.contour_paths(delta)
  
          # get vertices on boundary
          vertices = paths.vertices()
  
          # evaluate polynomial
          supremum = numpy.max(numpy.abs(p(vertices)))
  
          # compute bound
          bound.append(epsilon/(delta-epsilon)
                       * paths.length()/(2*numpy.pi*delta)
                       * supremum)
      return bound