docs: def interpolated_halfmax_points(x, y):
      """Given a curve y(x), find the x coordinates of points that have half the
      value of max(y), using linear interpolation. We're assuming that y(x) has
      a bandpass-ish shape, i.e., a single maximum and a drop to zero as we go
      to the edges of the function's domain. We also assume that x is sorted
      increasingly.
  
      """
      from scipy.interpolate import interp1d
      from scipy.optimize import fmin
  
      x = np.asarray(x)
      y = np.asarray(y)
      halfmax = 0.5 * y.max()
  
      # Guess from the actual samples.
  
      delta = y - halfmax
      guess1 = 0
      while delta[guess1] < 0:
          guess1 += 1
      guess2 = y.size - 1
      while delta[guess2] < 0:
          guess2 -= 1
  
      # Interpolate for fanciness.
  
      terp = interp1d(x, y, kind='linear', bounds_error=False, fill_value=0.)
      x1 = fmin(lambda x: (terp(x) - halfmax)**2, x[guess1], disp=False)
      x2 = fmin(lambda x: (terp(x) - halfmax)**2, x[guess2], disp=False)
  
      x1 = np.asscalar(x1)
      x2 = np.asscalar(x2)
  
      if x1 == x2:
          raise PKError('halfmax finding failed')
  
      if x1 > x2:
          x1, x2 = x2, x1
  
      return x1, x2