docs: def run_estimate_radius(data, dists, sample, d, rmin, rmax, ntry, run_parallel,
                          search_space='linspace', rad_bw_ratio=3.0):
      """
      This function is used to estimate the bandwidth, h, of the Gaussian Kernel:
          exp(-||x_i - x_j||/h^2)
  
  
      The "radius" refers to the truncation constant which we take to be
      rad_bw_ratio * h and this is the parameter over which we will iterate.
  
      We use the method of: https://arxiv.org/abs/1406.0118
  
      Parameters
      ----------
      data : numpy array,
          Original data set for which we are estimating the bandwidth
      dists : scipy.csr_matrix,
          A CSR matrix containing pairwise distances from data up to rmax
      sample : np.array,
          subset of data points over which to evaluate the radii
      d : int,
          dimension over which to evaluate the radii (smaller usually better)
      rmin : float,
          smallest radius ( = rad_bw_ratio * bandwidth) to consider
      rmax : float,
          largest radius ( = rad_bw_ratio * bandwidth) to consider
      ntry : int,
          number of radii between rmax and rmin to try
      run_parallel : bool,
          whether to run the analysis in parallel over radii
      search_space : str,
          either 'linspace' or 'logspace', choose to search in log or linear space
      rad_bw_ratio : float,
          the ratio of radius and kernel bandwidth, default to be 3 (radius = 3*h)
  
      Returns
      -------
      results : np.array (2d)
          first column is the set of radii and the second column is the distortion
          (smaller is better)
  
      """
      # declare global variables
      global DIST
      global X
  
      # process distance matrix
      dists = dists.tocsr()
      dists.data[dists.data > rmax] = 0.0
      dists.eliminate_zeros()
  
      # Assign global variables
      DIST = dists
      X = data
  
      n, dim = X.shape
      t0 = time.time()
      if run_parallel:
          ncpu = mp.cpu_count() # 32 for Newton
          processes = int(min(ntry, ncpu))
          print('using ' + str(processes) + ' processes to perform asynchronous parallel radius search')
          results = multi_process_radius_search(d, sample, rmin, rmax, ntry, processes, search_space, rad_bw_ratio)
      else:
          results = radius_search(d, sample, rmin, rmax, ntry, search_space, rad_bw_ratio)
      t1 = time.time()
      print('analysis took: ' + str(t1 - t0) + ' seconds to complete.')
      return(results)