docs: def RLS(anchors, W, r, print_out=False, grid=None, num_points=10):
      """ Range least squares (RLS) using grid search.
  
      Algorithm written by A.Beck, P.Stoica in "Approximate and Exact solutions of Source Localization Problems".
  
      :param anchors: anchor points
      :param r2: squared distances from anchors to point x.
      :param grid: where to search for solution.  (min, max) where min and max  are 
      lists of d elements, d being the dimension of the setup. If not given, the 
      search is conducted within the space covered by the anchors.
      :param num_points: number of grid points per direction.
  
      :return: estimated position of point x.
      """
      def cost_function(arr):
          X = np.c_[arr]
          r_measured = np.linalg.norm(anchors - X)
          mse = np.linalg.norm(r_measured - r)**2
          return mse
  
      if grid is None:
          grid = [np.min(anchors, axis=0), np.max(anchors, axis=0)]
  
      d = anchors.shape[1]
      x = np.linspace(grid[0][0], grid[1][0], num_points)
      y = np.linspace(grid[0][1], grid[1][1], num_points)
      if d == 2:
          errors_test = np.zeros((num_points, num_points))
          for i, xs in enumerate(x):
              for j, ys in enumerate(y):
                  errors_test[i, j] = cost_function((xs, ys))
          min_idx = errors_test.argmin()
          min_idx_multi = np.unravel_index(min_idx, errors_test.shape)
          xhat = np.c_[x[min_idx_multi[0]], y[min_idx_multi[1]]]
      elif d == 3:
          z = np.linspace(grid[0][2], grid[1][2], num_points)
          errors_test = np.zeros((num_points, num_points, num_points))
  
          # TODO: make this more efficient.
          #xx, yy, zz= np.meshgrid(x, y, z)
          for i, xs in enumerate(x):
              for j, ys in enumerate(y):
                  for k, zs in enumerate(z):
                      errors_test[i, j, k] = cost_function((xs, ys, zs))
          min_idx = errors_test.argmin()
          min_idx_multi = np.unravel_index(min_idx, errors_test.shape)
          xhat = np.c_[x[min_idx_multi[0]],
                       y[min_idx_multi[1]], z[min_idx_multi[2]]]
      else:
          raise ValueError('Non-supported number of dimensions.')
      return xhat[0]