docs: def compute_hulls(S, fS, domain):
      """
      (Re-)compute upper and lower hull given
      the segment points `S` with function values
      `fS` and the `domain` of the logpdf.
  
      Parameters
      ----------
      S : np.ndarray (N, 1)
         Straight-line segment points accumulated thus far.
  
      fS : tuple
          Value of the `logpdf` under sampling for each
          of the given segment points in `S`.
  
      domain : Tuple[float, float]
          Domain of `logpdf`.
          May be unbounded on either or both sides,
          in which case `(float("-inf"), float("inf"))`
          would be passed.
          If this domain is unbounded to the left,
          the derivative of the logpdf
          for x<= a must be positive.
          If this domain is unbounded to the right                  the derivative of the logpdf for x>=b
          must be negative.
  
      Returns
      ----------
      lower_hull: List[arspy.hull.HullNode]
      upper_hull: List[arspy.hull.HullNode]
  
      """
      assert(len(S) == len(fS))
      assert(len(domain) == 2)
  
      lower_hull = []
  
      for li in range(len(S) - 1):
          m = (fS[li + 1] - fS[li]) / (S[li + 1] - S[li])
          b = fS[li] - m * S[li]
          left = S[li]
          right = S[li + 1]
          lower_hull.append(HullNode(m=m, b=b, left=left, right=right))
  
      # compute upper piecewise-linear hull
  
      # expected final length of upper hull after full computation
      n_upper_segments = 2 * (len(S) - 2) + isinf(domain[0]) + isinf(domain[1])
  
      upper_hull = []
  
      if isinf(domain[0]):
          # first line (from -infinity)
  
          m = (fS[1] - fS[0]) / (S[1] - S[0])
          b = fS[0] - m * S[0]
          pr = compute_segment_log_prob(float("-inf"), S[0], m, b)
  
          upper_hull.append(
              HullNode(m=m, b=b, pr=pr, left=float("-inf"), right=S[0])
          )
  
      # second line
      m = (fS[2] - fS[1]) / (S[2] - S[1])
      b = fS[1] - m * S[1]
      pr = compute_segment_log_prob(S[0], S[1], m, b)
  
      upper_hull.append(HullNode(m=m, b=b, pr=pr, left=S[0], right=S[1]))
  
      # interior lines
      # there are two lines between each abscissa
      for li in range(1, len(S) - 2):
  
          m1 = (fS[li] - fS[li - 1]) / (S[li] - S[li - 1])
          b1 = fS[li] - m1 * S[li]
  
          m2 = (fS[li + 2] - fS[li + 1]) / (S[li + 2] - S[li + 1])
          b2 = fS[li + 1] - m2 * S[li + 1]
  
          if isinf(m1) and isinf(m2):
              raise ValueError("both hull slopes are infinite")
  
          dx1 = S[li] - S[li - 1]
          df1 = fS[li] - fS[li - 1]
          dx2 = S[li + 2] - S[li + 1]
          df2 = fS[li + 2] - fS[li + 1]
  
          f1 = fS[li]
          f2 = fS[li + 1]
          x1 = S[li]
          x2 = S[li + 1]
  
          # more numerically stable than above
          ix = ((f1 * dx1 - df1 * x1) * dx2 - (f2 * dx2 - df2 * x2) * dx1) / (df2 * dx1 - df1 * dx2)
  
          if isinf(m1) or abs(m1 - m2) < 10.0 ** 8 * eps(m1):
              ix = S[li]
              pr1 = float("-inf")
              pr2 = compute_segment_log_prob(ix, S[li + 1], m2, b2)
          elif isinf(m2):
              ix = S[li + 1]
              pr1 = compute_segment_log_prob(S[li], ix, m1, b1)
              pr2 = float("-inf")
          else:
              if isinf(ix):
                  raise ValueError("Non finite intersection")
  
              if abs(ix - S[li]) < 10.0 ** 12 * eps(S[li]):
                  ix = S[li]
              elif abs(ix - S[li + 1]) < 10.0**12 * eps(S[li + 1]):
                  ix = S[li + 1]
  
              if ix < S[li] or ix > S[li + 1]:
                  raise ValueError("Intersection out of bounds -- logpdf is not concave")
  
              pr1 = compute_segment_log_prob(S[li], ix, m1, b1)
              pr2 = compute_segment_log_prob(ix, S[li + 1], m2, b2)
  
              upper_hull.append(HullNode(m=m1, b=b1, pr=pr1, left=S[li], right=ix))
              upper_hull.append(HullNode(m=m2, b=b2, pr=pr2, left=ix, right=S[li + 1]))
  
      # second last line
      m = (fS[-2] - fS[-3]) / float(S[-2] - S[-3])
      b = fS[-2] - m * S[-2]
      pr = compute_segment_log_prob(S[-2], S[-1], m, b)
  
      upper_hull.append(HullNode(m=m, b=b, pr=pr, left=S[-2], right=S[-1]))
  
      if isinf(domain[1]):
          # last line (to infinity)
          m = (fS[-1] - fS[-2]) / (S[-1] - S[-2])
          b = fS[-1] - m * S[-1]
          pr = compute_segment_log_prob(S[-1], float("inf"), m, b)
  
          upper_hull.append(
              HullNode(m=m, b=b, pr=pr, left=S[-1], right=float("inf"))
          )
  
      # normalize probabilities
      normalized_probabilities = exp_normalize(
          asarray([node.pr for node in upper_hull])
      )
  
      for node, probability in zip(upper_hull, normalized_probabilities):
          node.pr = probability
  
      assert(len(lower_hull) == len(S) - 1)
      assert(len(upper_hull) == n_upper_segments)
  
      return lower_hull, upper_hull