code: def interpoled_resampling(W, x):
      
      N = W.shape[0] 
      idx = np.argsort(x)
      xs = x[idx] 
      ws = W[idx]
      cs = np.cumsum(avg_n_nplusone(ws)) 
      u = random.rand(N)
      xrs = np.empty(N)
      where = np.searchsorted(cs, u)
      # costs O(N log(N)) but algorithm has O(N log(N)) complexity anyway
      for n in range(N): 
          m = where[n]
          if m==0:
              xrs[n] = xs[0]
          elif m==N:
              xrs[n] = xs[-1]
          else:
              xrs[n] = interpol(cs[m-1], cs[m], xs[m-1], xs[m], u[n])
      return xrs
docs: Resampling based on an interpolated CDF, as described in Malik and Pitt. 
  
      Parameters
      ----------
      W: (N,) array
          weights 
      x: (N,) array 
          particles
  
      Returns
      -------
      xrs: (N,) array
          the resampled particles