docs: def dip_and_closest_unimodal_from_cdf(xF, yF, plotting=False, verbose=False, eps=1e-12):
      '''
          Dip computed as distance between empirical distribution function (EDF) and
          cumulative distribution function for the unimodal distribution with
          smallest such distance. The optimal unimodal distribution is found by
          the algorithm presented in
  
              Hartigan (1985): Computation of the dip statistic to test for
              unimodaliy. Applied Statistics, vol. 34, no. 3
  
          If the plotting option is enabled the optimal unimodal distribution
          function is plotted along with (xF, yF-dip) and (xF, yF+dip)
  
          xF  -   x-coordinates for EDF
          yF  -   y-coordinates for EDF
  
      '''
  
      ## TODO! Preprocess xF and yF so that yF increasing and xF does
      ## not have more than two copies of each x-value.
  
      if (xF[1:]-xF[:-1] < -eps).any():
          raise ValueError('Need sorted x-values to compute dip')
      if (yF[1:]-yF[:-1] < -eps).any():
          raise ValueError('Need sorted y-values to compute dip')
  
      # if plotting:
      #     Nplot = 5
      #     bfig = plt.figure(figsize=(12, 3))
      #     i = 1  # plot index
  
      D = 0  # lower bound for dip*2
  
      # [L, U] is interval where we still need to find unimodal function,
      # the modal interval
      L = 0
      U = len(xF) - 1
  
      # iGfin are the indices of xF where the optimal unimodal distribution is greatest
      # convex minorant to (xF, yF+dip)
      # iHfin are the indices of xF where the optimal unimodal distribution is least
      # concave majorant to (xF, yF-dip)
      iGfin = L
      iHfin = U
  
      while 1:
  
          iGG = greatest_convex_minorant_sorted(xF[L:(U+1)], yF[L:(U+1)])
          iHH = least_concave_majorant_sorted(xF[L:(U+1)], yF[L:(U+1)])
          iG = np.arange(L, U+1)[iGG]
          iH = np.arange(L, U+1)[iHH]
  
          # Interpolate. First and last point are in both and does not need
          # interpolation. Might cause trouble if included due to possiblity
          # of infinity slope at beginning or end of interval.
          if iG[0] != iH[0] or iG[-1] != iH[-1]:
              raise ValueError('Convex minorant and concave majorant should start and end at same points.')
          hipl = np.interp(xF[iG[1:-1]], xF[iH], yF[iH])
          gipl = np.interp(xF[iH[1:-1]], xF[iG], yF[iG])
          hipl = np.hstack([yF[iH[0]], hipl, yF[iH[-1]]])
          gipl = np.hstack([yF[iG[0]], gipl, yF[iG[-1]]])
          #hipl = lin_interpol_sorted(xF[iG], xF[iH], yF[iH])
          #gipl = lin_interpol_sorted(xF[iH], xF[iG], yF[iG])
  
          # Find largest difference between GCM and LCM.
          gdiff = hipl - yF[iG]
          hdiff = yF[iH] - gipl
          imaxdiffg = np.argmax(gdiff)
          imaxdiffh = np.argmax(hdiff)
          d = max(gdiff[imaxdiffg], hdiff[imaxdiffh])
  
          # # Plot current GCM and LCM.
          # if plotting:
          #     if i > Nplot:
          #         bfig = plt.figure(figsize=(12, 3))
          #         i = 1
          #     bax = bfig.add_subplot(1, Nplot, i)
          #     bax.plot(xF, yF, color='red')
          #     bax.plot(xF, yF-d/2, color='black')
          #     bax.plot(xF, yF+d/2, color='black')
          #     bax.plot(xF[iG], yF[iG]+d/2, color='blue')
          #     bax.plot(xF[iH], yF[iH]-d/2, color='blue')
  
          # if d <= D:
          #     if verbose:
          #         print("Difference in modal interval smaller than current dip")
          #     break
  
          # Find new modal interval so that largest difference is at endpoint
          # and set d to largest distance between current GCM and LCM.
          if gdiff[imaxdiffg] > hdiff[imaxdiffh]:
              L0 = iG[imaxdiffg]
              U0 = iH[iH >= L0][0]
          else:
              U0 = iH[imaxdiffh]
              L0 = iG[iG <= U0][-1]
          # Add points outside the modal interval to the final GCM and LCM.
          iGfin = np.hstack([iGfin, iG[(iG <= L0)*(iG > L)]])
          iHfin = np.hstack([iH[(iH >= U0)*(iH < U)], iHfin])
  
          # # Plot new modal interval
          # if plotting:
          #     ymin, ymax = bax.get_ylim()
          #     bax.axvline(xF[L0], ymin, ymax, color='orange')
          #     bax.axvline(xF[U0], ymin, ymax, color='red')
          #     bax.set_xlim(xF[L]-.1*(xF[U]-xF[L]), xF[U]+.1*(xF[U]-xF[L]))
  
          # Compute new lower bound for dip*2
          # i.e. largest difference outside modal interval
          gipl = np.interp(xF[L:(L0+1)], xF[iG], yF[iG])
          D = max(D, np.amax(yF[L:(L0+1)] - gipl))
          hipl = np.interp(xF[U0:(U+1)], xF[iH], yF[iH])
          D = max(D, np.amax(hipl - yF[U0:(U+1)]))
  
          if xF[U0]-xF[L0] < eps:
              if verbose:
                  print("Modal interval zero length")
              break
  
          # if plotting:
          #     mxpt = np.argmax(yF[L:(L0+1)] - gipl)
          #     bax.plot([xF[L:][mxpt], xF[L:][mxpt]], [yF[L:][mxpt]+d/2, 
          #       gipl[mxpt]+d/2], '+', color='red')
          #     mxpt = np.argmax(hipl - yF[U0:(U+1)])
          #     bax.plot([xF[U0:][mxpt], xF[U0:][mxpt]], [yF[U0:][mxpt]-d/2, 
          #       hipl[mxpt]-d/2], '+', color='red')
          #     i += 1
  
          # Change modal interval
          L = L0
          U = U0
  
          if d <= D:
              if verbose:
                  print("Difference in modal interval smaller than new dip")
              break
  
      # if plotting:
  
      #     # Add modal interval to figure
      #     bax.axvline(xF[L0], ymin, ymax, color='green', linestyle='dashed')
      #     bax.axvline(xF[U0], ymin, ymax, color='green', linestyle='dashed')
  
      #     ## Plot unimodal function (not distribution function)
      #     bfig = plt.figure()
      #     bax = bfig.add_subplot(1, 1, 1)
      #     bax.plot(xF, yF, color='red')
      #     bax.plot(xF, yF-D/2, color='black')
      #     bax.plot(xF, yF+D/2, color='black')
  
      # Find string position in modal interval
      iM = np.arange(iGfin[-1], iHfin[0]+1)
      yM_lower = yF[iM]-D/2
      yM_lower[0] = yF[iM[0]]+D/2
      iMM_concave = least_concave_majorant_sorted(xF[iM], yM_lower)
      iM_concave = iM[iMM_concave]
      #bax.plot(xF[iM], yM_lower, color='orange')
      #bax.plot(xF[iM_concave], yM_lower[iMM_concave], color='red')
      lcm_ipl = np.interp(xF[iM], xF[iM_concave], yM_lower[iMM_concave])
      try:
          mode = iM[np.nonzero(lcm_ipl > yF[iM]+D/2)[0][-1]]
          #bax.axvline(xF[mode], color='green', linestyle='dashed')
      except IndexError:
          iM_convex = np.zeros(0, dtype='i')
      else:
          after_mode = iM_concave > mode
          iM_concave = iM_concave[after_mode]
          iMM_concave = iMM_concave[after_mode]
          iM = iM[iM <= mode]
          iM_convex = iM[greatest_convex_minorant_sorted(xF[iM], yF[iM])]
  
      # if plotting:
  
      #     bax.plot(xF[np.hstack([iGfin, iM_convex, iM_concave, iHfin])],
      #              np.hstack([yF[iGfin] + D/2, yF[iM_convex] + D/2,
      #                         yM_lower[iMM_concave], yF[iHfin] - D/2]), color='blue')
      #     #bax.plot(xF[iM], yM_lower, color='orange')
  
      #     ## Plot unimodal distribution function
      #     bfig = plt.figure()
      #     bax = bfig.add_subplot(1, 1, 1)
      #     bax.plot(xF, yF, color='red')
      #     bax.plot(xF, yF-D/2, color='black')
      #     bax.plot(xF, yF+D/2, color='black')
  
      # Find string position in modal interval
      iM = np.arange(iGfin[-1], iHfin[0]+1)
      yM_lower = yF[iM]-D/2
      yM_lower[0] = yF[iM[0]]+D/2
      iMM_concave = least_concave_majorant_sorted(xF[iM], yM_lower)
      iM_concave = iM[iMM_concave]
      #bax.plot(xF[iM], yM_lower, color='orange')
      #bax.plot(xF[iM_concave], yM_lower[iMM_concave], color='red')
      lcm_ipl = np.interp(xF[iM], xF[iM_concave], yM_lower[iMM_concave])
      try:
          mode = iM[np.nonzero(lcm_ipl > yF[iM]+D/2)[0][-1]]
          #bax.axvline(xF[mode], color='green', linestyle='dashed')
      except IndexError:
          iM_convex = np.zeros(0, dtype='i')
      else:
          after_mode = iM_concave > mode
          iM_concave = iM_concave[after_mode]
          iMM_concave = iMM_concave[after_mode]
          iM = iM[iM <= mode]
          iM_convex = iM[greatest_convex_minorant_sorted(xF[iM], yF[iM])]
  
      # Closest unimodal curve
      xU = xF[np.hstack([iGfin[:-1], iM_convex, iM_concave, iHfin[1:]])]
      yU = np.hstack([yF[iGfin[:-1]] + D/2, yF[iM_convex] + D/2,
                      yM_lower[iMM_concave], yF[iHfin[1:]] - D/2])
      # Add points so unimodal curve goes from 0 to 1
      k_start = (yU[1]-yU[0])/(xU[1]-xU[0]+1e-5)
      xU_start = xU[0] - yU[0]/(k_start+1e-5)
      k_end = (yU[-1]-yU[-2])/(xU[-1]-xU[-2]+1e-5)
      xU_end = xU[-1] + (1-yU[-1])/(k_end+1e-5)
      xU = np.hstack([xU_start, xU, xU_end])
      yU = np.hstack([0, yU, 1])
  
      # if plotting:
      #     bax.plot(xU, yU, color='blue')
      #     #bax.plot(xF[iM], yM_lower, color='orange')
      #     plt.show()
  
      return D/2, (xU, yU)