docs: def S2_surface(self, sizes, bounds, presets, covers, use_torch=False,
              num_samples = 10):
          """Calculates the sensitivity surface of a GrFN for the two variables with
          the highest S2 index.
  
          Args:
              num_samples: Number of samples for sensitivity analysis.
              sizes: Tuple of (number of x inputs, number of y inputs).
              bounds: Set of bounds for GrFN inputs.
              presets: Set of standard values for GrFN inputs.
  
          Returns:
              Tuple:
                  Tuple: The names of the two variables that were selected
                  Tuple: The X, Y vectors of eval values
                  Z: The numpy matrix of output evaluations
  
          """
          args = self.inputs
          Si = self.sobol_analysis(
              num_samples,
              {
                  "num_vars": len(args),
                  "names": args,
                  "bounds": [bounds[arg] for arg in args],
              },
              covers
          )
          S2 = Si["S2"]
          (s2_max, v1, v2) = get_max_s2_sensitivity(S2)
  
          x_var = args[v1]
          y_var = args[v2]
          search_space = [(x_var, bounds[x_var]), (y_var, bounds[y_var])]
          preset_vals = {
              arg: presets[arg]
              for i, arg in enumerate(args)
              if i != v1 and i != v2
          }
  
          X = np.linspace(*search_space[0][1], sizes[0])
          Y = np.linspace(*search_space[1][1], sizes[1])
  
          if use_torch:
              Xm, Ym = torch.meshgrid(torch.tensor(X), torch.tensor(Y))
              inputs = {n: torch.full_like(Xm, v) for n, v in presets.items()}
              inputs.update({search_space[0][0]: Xm, search_space[1][0]: Ym})
              Z = self.run(inputs, covers).numpy()
          else:
              Xm, Ym = np.meshgrid(X, Y)
              Z = np.zeros((len(X), len(Y)))
              for x, y in itertools.product(range(len(X)), range(len(Y))):
                  inputs = {n: v for n, v in presets.items()}
                  inputs.update({search_space[0][0]: x, search_space[1][0]: y})
                  Z[x][y] = self.run(inputs, covers)
  
          return X, Y, Z, x_var, y_var