code: def partial_dependence_2D(space, model, i, j, samples, n_points=40):
      """
      Calculate the partial dependence for two dimensions in the search-space.
  
      This uses the given model to calculate the average objective value
      for all the samples, where the given dimensions are fixed at
      regular intervals between their bounds.
  
      This shows how the given dimensions affect the objective value
      when the influence of all other dimensions are averaged out.
  
      Parameters
      ----------
      space : `Space`
          The parameter space over which the minimization was performed.
  
      model
          Surrogate model for the objective function.
  
      i : int
          The first dimension for which to calculate the partial dependence.
  
      j : int
          The second dimension for which to calculate the partial dependence.
  
      samples : np.array, shape=(n_points, n_dims)
          Randomly sampled and transformed points to use when averaging
          the model function at each of the `n_points` when using partial
          dependence.
  
      n_points : int, default=40
          Number of points at which to evaluate the partial dependence
          along each dimension `i` and `j`.
  
      Returns
      -------
      xi : np.array, shape=n_points
          The points at which the partial dependence was evaluated.
  
      yi : np.array, shape=n_points
          The points at which the partial dependence was evaluated.
  
      zi : np.array, shape=(n_points, n_points)
          The average value of the objective function at each point `(xi, yi)`.
      """
      # The idea is to step through one dimension, evaluating the model with
      # that dimension fixed and averaging either over random values or over
      # the given ones in x_val in all other dimensions.
      # (Or step through 2 dimensions when i and j are given.)
      # Categorical dimensions make this interesting, because they are one-
      # hot-encoded, so there is a one-to-many mapping of input dimensions
      # to transformed (model) dimensions.
  
      # dim_locs[i] is the (column index of the) start of dim i in
      # sample_points.
      # This is usefull when we are using one hot encoding, i.e using
      # categorical values
      dim_locs = np.cumsum(
          [0] + [d.transformed_size for d in space.dimensions]
      )
  
      def _calc(x, y):
          """
          Helper-function to calculate the average predicted
          objective value for the given model, when setting
          the index1'th dimension of the search-space to the value x
          and setting the index2'th dimension to the value y,
          and then averaging over all samples.
          """
          rvs_ = np.array(samples)  # copy
          rvs_[:, dim_locs[j] : dim_locs[j + 1]] = x
          rvs_[:, dim_locs[i] : dim_locs[i + 1]] = y
          return np.mean(model.predict(rvs_))
  
      xi, xi_transformed = _evenly_sample(space.dimensions[j], n_points)
      yi, yi_transformed = _evenly_sample(space.dimensions[i], n_points)
      # Calculate the partial dependence for all combinations of these points.
      zi = [
          [_calc(x, y) for x in xi_transformed] for y in yi_transformed
      ]
  
      # Convert list-of-list to a numpy array.
      zi = np.array(zi)
  
      return xi, yi, zi
  
imports: import sys
  import numpy
  from itertools import count
  from functools import partial
  from scipy.optimize import OptimizeResult
  from .acquisition import _gaussian_acquisition
  from guild._skopt import expected_minimum, expected_minimum_random_sampling
  from .space import Categorical
  from collections import Counter
  import matplotlib.pyplot
  from matplotlib.pyplot import cm
  from matplotlib.ticker import LogLocator
  from matplotlib.ticker import MaxNLocator, FuncFormatter
label_desc: Use image node_id 10 for calling a global function with example usage: partial_dependence_2D(space, model, i, j, samples, n_points) and returns: xi,  yi,  zi, np