Skip to main content
`Ctrl`+`K`

# Source code for skopt.plots

    """Plotting functions."""

    import sys
    import warnings
    from collections import Counter
    from collections.abc import Iterable
    from functools import partial
    from itertools import count

    import numpy as np
    from scipy.optimize import OptimizeResult

    from skopt import expected_minimum, expected_minimum_random_sampling

    from .acquisition import _gaussian_acquisition
    from .space import Categorical

    # For plot tests, matplotlib must be set to headless mode early
    if 'pytest' in sys.modules:
        import matplotlib

        matplotlib.use('Agg')

    from matplotlib import gridspec
    from matplotlib import pyplot as plt
    from matplotlib.pyplot import cm
    from matplotlib.ticker import FuncFormatter, LogLocator, MaxNLocator  # noqa: E402

    def plot_convergence(*args, true_minimum=None, yscale=None, ax=None):
        """Plot one or several convergence traces.

        Parameters
        ----------
        results: `OptimizeResult`, iterable of `OptimizeResult`, or a 2-tuple
            of a label and a `OptimizeResult` or an iterable of `OptimizeResult`.
            The result(s) for which to plot the convergence trace.

            - if an `OptimizeResult`, draw the corresponding single trace
            - if an iterable of `OptimizeResult`, draw all traces in the same
              plot as well as the average convergence trace
            - if a tuple, the label names the trace(s) and the behavior is as
              specified above.

        true_minimum : float, optional
            The true minimum value of the function, if known.

        yscale : None or string, optional
            The scale for the y-axis.

        Returns
        -------
        ax : `Axes`
            The matplotlib axes the plot was drawn in
        """

        if ax is None:
            ax = plt.gca()

        ax.set_title("Convergence plot")
        ax.set_xlabel("Number of calls $n$")
        ax.set_ylabel(r"$\min f(x)$ after $n$ calls")
        ax.grid()

        if yscale is not None:
            ax.set_yscale(yscale)

        colors = cm.viridis(np.linspace(0.25, 1.0, len(args)))

        for index, (arg, color) in enumerate(zip(args, colors)):
            if isinstance(arg, tuple):
                label, arg = arg
            else:
                label = None

            if isinstance(arg, OptimizeResult):
                opt_res = arg
                n_calls = len(opt_res.x_iters)
                mins = [np.min(opt_res.func_vals[: i + 1]) for i in range(n_calls)]
                ax.plot(
                    range(1, n_calls + 1),
                    mins,
                    c=color,
                    marker=".",
                    markersize=12,
                    lw=2,
                    label=label,
                )

            elif isinstance(arg, Iterable) and all(
                isinstance(elem, OptimizeResult) for elem in arg
            ):
                mins = [
                    [
                        np.min(opt_res.func_vals[: i + 1])
                        for i in range(len(opt_res.x_iters))
                    ]
                    for opt_res in arg
                ]

                # graciously handle "ragged" array, i.e. differing length arrays
                max_n_calls = max(len(m) for m in mins)
                mean_arr = np.empty((len(mins), max_n_calls))
                mean_arr[:] = np.nan

                for i, m in enumerate(mins):
                    ax.plot(range(1, 1 + len(m)), m, c=color, alpha=0.2)
                    mean_arr[i, : len(m)] = m

                if np.isnan(mean_arr).any():
                    warnings.warn(
                        "Inconsistent number of function calls in "
                        f"argument at pos {index}"
                    )

                ax.plot(
                    range(1, 1 + max_n_calls),
                    np.nanmean(mins, axis=0),
                    c=color,
                    marker=".",
                    markersize=12,
                    lw=2,
                    label=label,
                )

            else:
                raise ValueError(
                    "Cannot plot convergence trace for "
                    f"{arg.__class__.__name__} object {arg}"
                )

        if true_minimum:
            ax.axhline(true_minimum, linestyle="--", color="r", lw=1, label="True minimum")

        if true_minimum or label:
            ax.legend(loc="best")

        return ax

    def plot_gaussian_process(
        res,
        ax=None,
        n_calls=-1,
        objective=None,
        n_points=1000,
        noise_level=0,
        show_legend=True,
        show_title=True,
        show_acq_func=False,
        show_next_point=False,
        show_observations=True,
        show_mu=True,
    ):
        """Plots the optimization results and the gaussian process for 1-D objective
        functions.

        Parameters
        ----------
        res :  `OptimizeResult`
            The result for which to plot the gaussian process.

        ax : `Axes`, optional
            The matplotlib axes on which to draw the plot, or `None` to create
            a new one.

        n_calls : int, default: -1
            Can be used to evaluate the model at call `n_calls`.

        objective : func, default: None
            Defines the true objective function. Must have one input parameter.

        n_points : int, default: 1000
            Number of data points used to create the plots

        noise_level : float, default: 0
            Sets the estimated noise level

        show_legend : boolean, default: True
            When True, a legend is plotted.

        show_title : boolean, default: True
            When True, a title containing the found minimum value
            is shown

        show_acq_func : boolean, default: False
            When True, the acquisition function is plotted

        show_next_point : boolean, default: False
            When True, the next evaluated point is plotted

        show_observations : boolean, default: True
            When True, observations are plotted as dots.

        show_mu : boolean, default: True
            When True, the predicted model is shown.

        Returns
        -------
        ax : `Axes`
            The matplotlib axes.
        """

        if ax is None:
            ax = plt.gca()
        n_dims = res.space.n_dims
        assert n_dims == 1, "Space dimension must be 1"
        dimension = res.space.dimensions[0]
        x, x_model = _evenly_sample(dimension, n_points)
        x = x.reshape(-1, 1)
        x_model = x_model.reshape(-1, 1)
        if res.specs is not None and "args" in res.specs:
            n_random = res.specs["args"].get('n_random_starts', None)
            acq_func = res.specs["args"].get("acq_func", "EI")
            acq_func_kwargs = res.specs["args"].get("acq_func_kwargs", {})

        if acq_func_kwargs is None:
            acq_func_kwargs = {}
        if acq_func is None or acq_func == "gp_hedge":
            acq_func = "EI"
        if n_random is None:
            n_random = len(res.x_iters) - len(res.models)

        if objective is not None:
            fx = np.array([objective(x_i) for x_i in x])
        if n_calls < 0:
            model = res.models[-1]
            curr_x_iters = res.x_iters
            curr_func_vals = res.func_vals
        else:
            model = res.models[n_calls]

            curr_x_iters = res.x_iters[: n_random + n_calls]
            curr_func_vals = res.func_vals[: n_random + n_calls]

        # Plot true function.
        if objective is not None:
            ax.plot(x, fx, "r--", label="True (unknown)")
            ax.fill(
                np.concatenate([x, x[::-1]]),
                np.concatenate(
                    (
                        [fx_i - 1.9600 * noise_level for fx_i in fx],
                        [fx_i + 1.9600 * noise_level for fx_i in fx[::-1]],
                    )
                ),
                alpha=0.2,
                fc="r",
                ec="None",
            )

        # Plot GP(x) + contours
        if show_mu:
            per_second = acq_func.endswith("ps")
            if per_second:
                y_pred, sigma = model.estimators_[0].predict(x_model, return_std=True)
            else:
                y_pred, sigma = model.predict(x_model, return_std=True)
            ax.plot(x, y_pred, "g--", label=r"$\mu_{GP}(x)$")
            ax.fill(
                np.concatenate([x, x[::-1]]),
                np.concatenate([y_pred - 1.9600 * sigma, (y_pred + 1.9600 * sigma)[::-1]]),
                alpha=0.2,
                fc="g",
                ec="None",
            )

        # Plot sampled points
        if show_observations:
            ax.plot(curr_x_iters, curr_func_vals, "r.", markersize=8, label="Observations")
        if (show_mu or show_observations or objective is not None) and show_acq_func:
            ax_ei = ax.twinx()
            ax_ei.set_ylabel(str(acq_func) + "(x)")
            plot_both = True
        else:
            ax_ei = ax
            plot_both = False
        if show_acq_func:
            acq = _gaussian_acquisition(
                x_model,
                model,
                y_opt=np.min(curr_func_vals),
                acq_func=acq_func,
                acq_func_kwargs=acq_func_kwargs,
            )
            next_x = x[np.argmin(acq)]
            next_acq = acq[np.argmin(acq)]
            acq = -acq
            next_acq = -next_acq
            ax_ei.plot(x, acq, "b", label=str(acq_func) + "(x)")
            if not plot_both:
                ax_ei.fill_between(x.ravel(), 0, acq.ravel(), alpha=0.3, color='blue')

            if show_next_point and next_x is not None:
                ax_ei.plot(next_x, next_acq, "bo", markersize=6, label="Next query point")

        if show_title:
            ax.set_title(fr"x* = {res.x[0]:.4f}, f(x*) = {res.fun:.4f}")
        # Adjust plot layout
        ax.grid()
        ax.set_xlabel("x")
        ax.set_ylabel("f(x)")
        if show_legend:
            if plot_both:
                lines, labels = ax.get_legend_handles_labels()
                lines2, labels2 = ax_ei.get_legend_handles_labels()
                ax_ei.legend(
                    lines + lines2,
                    labels + labels2,
                    loc="best",
                    prop={'size': 6},
                    numpoints=1,
                )
            else:
                ax.legend(loc="best", prop={'size': 6}, numpoints=1)

        return ax

    def plot_regret(*args, ax=None, true_minimum=None, yscale=None):
        """Plot one or several cumulative regret traces.

        Parameters
        ----------
        args[i] : `OptimizeResult`, list of `OptimizeResult`, or tuple
            The result(s) for which to plot the cumulative regret trace.

            - if `OptimizeResult`, then draw the corresponding single trace;
            - if list of `OptimizeResult`, then draw the corresponding cumulative
                regret traces in transparency, along with the average cumulative
                regret trace;
            - if tuple, then `args[i][0]` should be a string label and `args[i][1]`
              an `OptimizeResult` or a list of `OptimizeResult`.

        ax : Axes`, optional
            The matplotlib axes on which to draw the plot, or `None` to create
            a new one.

        true_minimum : float, optional
            The true minimum value of the function, if known.

        yscale : None or string, optional
            The scale for the y-axis.

        Returns
        -------
        ax : `Axes`
            The matplotlib axes.
        """

        if ax is None:
            ax = plt.gca()

        ax.set_title("Cumulative regret plot")
        ax.set_xlabel("Number of calls $n$")
        ax.set_ylabel(r"$\sum_{i=0}^n(f(x_i) - optimum)$ after $n$ calls")
        ax.grid()

        if yscale is not None:
            ax.set_yscale(yscale)

        colors = cm.viridis(np.linspace(0.25, 1.0, len(args)))

        if true_minimum is None:
            results = []
            for res in args:
                if isinstance(res, tuple):
                    res = res[1]

                if isinstance(res, OptimizeResult):
                    results.append(res)
                elif isinstance(res, list):
                    results.extend(res)
            true_minimum = np.min([np.min(r.func_vals) for r in results])

        for results, color in zip(args, colors):
            if isinstance(results, tuple):
                name, results = results
            else:
                name = None

            if isinstance(results, OptimizeResult):
                n_calls = len(results.x_iters)
                regrets = [
                    np.sum(results.func_vals[:i] - true_minimum)
                    for i in range(1, n_calls + 1)
                ]
                ax.plot(
                    range(1, n_calls + 1),
                    regrets,
                    c=color,
                    marker=".",
                    markersize=12,
                    lw=2,
                    label=name,
                )

            elif isinstance(results, list):
                n_calls = len(results[0].x_iters)
                iterations = range(1, n_calls + 1)
                regrets = [
                    [np.sum(r.func_vals[:i] - true_minimum) for i in iterations]
                    for r in results
                ]

                for cr in regrets:
                    ax.plot(iterations, cr, c=color, alpha=0.2)

                ax.plot(
                    iterations,
                    np.mean(regrets, axis=0),
                    c=color,
                    marker=".",
                    markersize=12,
                    lw=2,
                    label=name,
                )

        if name:
            ax.legend(loc="best")

        return ax

    def _format_scatter_plot_axes(ax, space, ylabel, plot_dims, dim_labels=None):
        # Work out min, max of y axis for the diagonal so we can adjust
        # them all to the same value
        diagonal_ylim = _get_ylim_diagonal(ax)

        # Number of search-space dimensions we are using.
        if isinstance(ax, (list, np.ndarray)):
            n_dims = len(plot_dims)
        else:
            n_dims = 1

        if dim_labels is None:
            dim_labels = [
                "$X_{%i}$" % i if d.name is None else d.name for i, d in plot_dims
            ]
        # Axes for categorical dimensions are really integers; we have to
        # label them with the category names
        iscat = [isinstance(dim[1], Categorical) for dim in plot_dims]

        # Deal with formatting of the axes
        for i in range(n_dims):  # rows
            for j in range(n_dims):  # columns
                if isinstance(ax, np.ndarray):
                    ax_ = ax[i, j]
                else:
                    ax_ = ax
                index_i, dim_i = plot_dims[i]
                index_j, dim_j = plot_dims[j]
                if j > i:
                    ax_.axis("off")
                elif i > j:  # off-diagonal plots
                    # plots on the diagonal are special, like Texas. They have
                    # their own range so do not mess with them.
                    if not iscat[i]:  # bounds not meaningful for categoricals
                        ax_.set_ylim(*dim_i.bounds)
                    if iscat[j]:
                        # partial() avoids creating closures in a loop
                        ax_.xaxis.set_major_formatter(
                            FuncFormatter(partial(_cat_format, dim_j))
                        )
                    else:
                        ax_.set_xlim(*dim_j.bounds)
                    if j == 0:  # only leftmost column (0) gets y labels
                        ax_.set_ylabel(dim_labels[i])
                        if iscat[i]:  # Set category labels for left column
                            ax_.yaxis.set_major_formatter(
                                FuncFormatter(partial(_cat_format, dim_i))
                            )
                    else:
                        ax_.set_yticklabels([])

                    # for all rows except ...
                    if i < n_dims - 1:
                        ax_.set_xticklabels([])
                    # ... the bottom row
                    else:
                        [label.set_rotation(45) for label in ax_.get_xticklabels()]
                        ax_.set_xlabel(dim_labels[j])

                    # configure plot for linear vs log-scale
                    if dim_j.prior == 'log-uniform':
                        ax_.set_xscale('log')
                    else:
                        ax_.xaxis.set_major_locator(
                            MaxNLocator(6, prune='both', integer=iscat[j])
                        )

                    if dim_i.prior == 'log-uniform':
                        ax_.set_yscale('log')
                    else:
                        ax_.yaxis.set_major_locator(
                            MaxNLocator(6, prune='both', integer=iscat[i])
                        )

                else:  # diagonal plots
                    ax_.set_ylim(*diagonal_ylim)
                    if not iscat[i]:
                        low, high = dim_i.bounds
                        ax_.set_xlim(low, high)
                    ax_.yaxis.tick_right()
                    ax_.yaxis.set_label_position('right')
                    ax_.yaxis.set_ticks_position('both')
                    ax_.set_ylabel(ylabel)

                    ax_.xaxis.tick_top()
                    ax_.xaxis.set_label_position('top')
                    ax_.set_xlabel(dim_labels[j])

                    if dim_i.prior == 'log-uniform':
                        ax_.set_xscale('log')
                    else:
                        ax_.xaxis.set_major_locator(
                            MaxNLocator(6, prune='both', integer=iscat[i])
                        )
                        if iscat[i]:
                            ax_.xaxis.set_major_formatter(
                                FuncFormatter(partial(_cat_format, dim_i))
                            )

        return ax

    def _make_subgrid(ax, n_rows, n_cols=None, fig_kwargs_=None, **gridspec_kwargs):
        """Makes a subgrid inside an existing axis object."""
        if n_cols is None:
            n_cols = n_rows
        fig_kwargs_ = fig_kwargs_ or {}
        if ax is None:
            fig, ax = plt.subplots(**fig_kwargs_)
        else:
            fig = ax.get_figure()

        grid_spec = gridspec.GridSpecFromSubplotSpec(
            n_rows, n_cols, subplot_spec=ax.get_subplotspec(), **gridspec_kwargs  # noqa
        )
        axes = np.empty((n_rows, n_cols), dtype=object)
        for i in range(n_rows):
            for j in range(n_cols):
                axes[i, j] = fig.add_subplot(grid_spec[i, j])
        ax.get_xaxis().set_visible(False)
        ax.get_yaxis().set_visible(False)
        for spine in ax.spines.values():
            spine.set_visible(False)
        return ax, axes

    def partial_dependence(
        space, model, i, j=None, sample_points=None, n_samples=250, n_points=40, x_eval=None
    ):
        """Calculate the partial dependence for dimensions `i` and `j` with respect to the
        objective value, as approximated by `model`.

        The partial dependence plot shows how the value of the dimensions
        `i` and `j` influence the `model` predictions after "averaging out"
        the influence of all other dimensions.

        When `x_eval` is not `None`, the given values are used instead of
        random samples. In this case, `n_samples` will be ignored.

        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, default=None
            The second dimension for which to calculate the partial dependence.
            To calculate the 1D partial dependence on `i` alone set `j=None`.

        sample_points : np.array, shape=(n_points, n_dims), default=None
            Only used when `x_eval=None`, i.e in case partial dependence should
            be calculated.
            Randomly sampled and transformed points to use when averaging
            the model function at each of the `n_points` when using partial
            dependence.

        n_samples : int, default=100
            Number of random samples to use for averaging the model function
            at each of the `n_points` when using partial dependence. Only used
            when `sample_points=None` and `x_eval=None`.

        n_points : int, default=40
            Number of points at which to evaluate the partial dependence
            along each dimension `i` and `j`.

        x_eval : list, default=None
            `x_eval` is a list of parameter values or None. In case `x_eval`
            is not None, the parsed dependence will be calculated using these
            values.
            Otherwise, random selected samples will be used.

        Returns
        -------
        For 1D partial dependence:

        xi : np.array
            The points at which the partial dependence was evaluated.

        yi : np.array
            The value of the model at each point `xi`.

        For 2D partial dependence:

        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 value of the model at each point `(xi, yi)`.

        For Categorical variables, the `xi` (and `yi` for 2D) returned are
        the indices of the variable in `Dimension.categories`.
        """
        # If we haven't parsed an x_eval list we use random sampled values instead
        if x_eval is None and sample_points is None:
            sample_points = space.transform(space.rvs(n_samples=n_samples))
        elif sample_points is None:
            sample_points = space.transform([x_eval])

        if j is None:
            return partial_dependence_1D(space, model, i, sample_points, n_points)
        else:
            return partial_dependence_2D(space, model, i, j, sample_points, n_points)

    def plot_objective(
        result,
        levels=10,
        n_points=40,
        n_samples=250,
        size=2,
        wspace=0.35,
        hspace=0.35,
        zscale='linear',
        dimensions=None,
        sample_source='random',
        minimum='result',
        n_minimum_search=None,
        plot_dims=None,
        show_points=True,
        cmap='viridis_r',
        ax=None,
    ):
        """Plot a 2-d matrix with so-called Partial Dependence plots of the objective
        function. This shows the influence of each search-space dimension on the objective
        function.

        This uses the last fitted model for estimating the objective function.

        The diagonal shows the effect of a single dimension on the
        objective function, while the plots below the diagonal show
        the effect on the objective function when varying two dimensions.

        The Partial Dependence is calculated by averaging the objective value
        for a number of random samples in the search-space,
        while keeping one or two dimensions fixed at regular intervals. This
        averages out the effect of varying the other dimensions and shows
        the influence of one or two dimensions on the objective function.

        Also shown are small black dots for the points that were sampled
        during optimization.

        A red star indicates per default the best observed minimum, but
        this can be changed by changing argument ´minimum´.

        .. note::
              The Partial Dependence plot is only an estimation of the surrogate
              model which in turn is only an estimation of the true objective
              function that has been optimized. This means the plots show
              an "estimate of an estimate" and may therefore be quite imprecise,
              especially if few samples have been collected during the
              optimization
              (e.g. less than 100-200 samples), and in regions of the search-space
              that have been sparsely sampled (e.g. regions away from the optimum).
              This means that the plots may change each time you run the
              optimization and they should not be considered completely reliable.
              These compromises are necessary because we cannot evaluate the
              expensive objective function in order to plot it, so we have to use
              the cheaper surrogate model to plot its contour. And in order to
              show search-spaces with 3 dimensions or more in a 2-dimensional
              plot,
              we further need to map those dimensions to only 2-dimensions using
              the Partial Dependence, which also causes distortions in the plots.

        Parameters
        ----------
        result : `OptimizeResult`
            The optimization results from calling e.g. `gp_minimize()`.

        levels : int, default=10
            Number of levels to draw on the contour plot, passed directly
            to `plt.contourf()`.

        n_points : int, default=40
            Number of points at which to evaluate the partial dependence
            along each dimension.

        n_samples : int, default=250
            Number of samples to use for averaging the model function
            at each of the `n_points` when `sample_method` is set to 'random'.

        size : float, default=2
            Height (in inches) of each facet. Ignored if ``ax`` is provided.

        wspace : float, default=0.35
            The width of the padding between subplots, as a fraction of the
            average Axes width. Ignored if ``ax`` is provided.

        hspace : float, default=0.35
           The height of the padding between subplots, as a fraction of the
           average Axes height. Ignored if ``ax`` is provided.

        zscale : str, default='linear'
            Scale to use for the z axis of the contour plots. Either 'linear'
            or 'log'.

        dimensions : list of str, default=None
            Labels of the dimension
            variables. `None` defaults to `space.dimensions[i].name`, or
            if also `None` to `['X_0', 'X_1', ..]`.

        plot_dims : list of str and int, default=None
            List of dimension names or dimension indices from the
            search-space dimensions to be included in the plot.
            If `None` then use all dimensions except constant ones
            from the search-space.

        sample_source : str or list of floats, default='random'
            Defines to samples generation to use for averaging the model function
            at each of the `n_points`.

            A partial dependence plot is only generated, when `sample_source`
            is set to 'random' and `n_samples` is sufficient.

            `sample_source` can also be a list of
            floats, which is then used for averaging.

            Valid strings:

            - 'random' - `n_samples` random samples will used
            - 'result' - Use only the best observed parameters
            - 'expected_minimum' - Parameters that gives the best
              minimum Calculated using scipy's minimize method.
              This method currently does not work with categorical values.
            - 'expected_minimum_random' - Parameters that gives the
              best minimum when using naive random sampling.
              Works with categorical values.

        minimum : str or list of floats, default = 'result'
            Defines the values for the red points in the plots.
            Valid strings:

            - 'result' - Use best observed parameters
            - 'expected_minimum' - Parameters that gives the best
              minimum Calculated using scipy's minimize method.
              This method currently does not work with categorical values.
            - 'expected_minimum_random' - Parameters that gives the
              best minimum when using naive random sampling.
              Works with categorical values

        n_minimum_search : int, default = None
            Determines how many points should be evaluated
            to find the minimum when using 'expected_minimum' or
            'expected_minimum_random'. Parameter is used when
            `sample_source` and/or `minimum` is set to
            'expected_minimum' or 'expected_minimum_random'.

        show_points: bool, default = True
            Choose whether to show evaluated points in the
            contour plots.

        cmap: str or Colormap, default = 'viridis_r'
            Color map for contour plots. Passed directly to
            `plt.contourf()`

        ax: `Matplotlib.Axes`, default= None
            An axis object in which to plot the dependence plot. If provided,
            ``size`` is ignored and the caller is responsible for the size of the
            plot.

        Returns
        -------
        ax : `Matplotlib.Axes`
            The axes object the plot was drawn in
        """
        # Here we define the values for which to plot the red dot (2d plot) and
        # the red dotted line (1d plot).
        # These same values will be used for evaluating the plots when
        # calculating dependence. (Unless partial
        # dependence is to be used instead).
        space = result.space
        # Get the relevant search-space dimensions.
        if plot_dims is None:
            # Get all dimensions.
            plot_dims = []
            for row in range(space.n_dims):
                if space.dimensions[row].is_constant:
                    continue
                plot_dims.append((row, space.dimensions[row]))
        else:
            plot_dims = space[plot_dims]
        # Number of search-space dimensions we are using.
        n_dims = len(plot_dims)
        if dimensions is not None:
            assert len(dimensions) == n_dims
        x_vals = _evaluate_min_params(result, minimum, n_minimum_search)
        if sample_source == "random":
            x_eval = None
            samples = space.transform(space.rvs(n_samples=n_samples))
        else:
            x_eval = _evaluate_min_params(result, sample_source, n_minimum_search)
            samples = space.transform([x_eval])
        x_samples, minimum, _ = _map_categories(space, result.x_iters, x_vals)

        if zscale == 'log':
            locator = LogLocator()
        elif zscale == 'linear':
            locator = None
        else:
            raise ValueError(
                "Valid values for zscale are 'linear' and 'log'," " not '%s'." % zscale
            )

        fig_kwargs = dict(figsize=(size * n_dims, size * n_dims))
        ax, axes = _make_subgrid(
            ax, n_dims, fig_kwargs_=fig_kwargs, wspace=wspace, hspace=hspace
        )

        for i in range(n_dims):
            for j in range(n_dims):
                if i == j:
                    index, _ = plot_dims[i]
                    xi, yi = partial_dependence_1D(
                        space, result.models[-1], index, samples=samples, n_points=n_points
                    )
                    ax_ = axes[i, j]
                    ax_.plot(xi, yi)
                    ax_.axvline(minimum[index], linestyle="--", color="r", lw=1)

                # lower triangle
                elif i > j:
                    index1, _ = plot_dims[i]
                    index2, _ = plot_dims[j]
                    xi, yi, zi = partial_dependence_2D(
                        space, result.models[-1], index1, index2, samples, n_points
                    )
                    ax_ = axes[i, j]
                    ax_.contourf(xi, yi, zi, levels, locator=locator, cmap=cmap)
                    if show_points:
                        ax_.scatter(
                            x_samples[:, index2], x_samples[:, index1], c='k', s=10, lw=0.0
                        )
                    ax_.scatter(
                        minimum[index2], minimum[index1], c=['r'], s=100, lw=0.0, marker='*'
                    )
        ylabel = "Partial dependence"

        # Make various adjustments to the plots.
        _format_scatter_plot_axes(
            axes, space, ylabel=ylabel, plot_dims=plot_dims, dim_labels=dimensions
        )
        return ax

    def plot_evaluations(
        result,
        bins=20,
        dimensions=None,
        plot_dims=None,
        size=2,
        wspace=0.35,
        hspace=0.35,
        cmap="viridis",
        ax=None,
    ):
        """Visualize the order in which points were sampled during optimization.

        This creates a 2-d matrix plot where the diagonal plots are histograms
        that show the distribution of samples for each search-space dimension.

        The plots below the diagonal are scatter-plots of the samples for
        all combinations of search-space dimensions.

        The order in which samples
        were evaluated is encoded in each point's color.

        A red star shows the best found parameters.

        Parameters
        ----------
        result : `OptimizeResult`
            The optimization results from calling e.g. `gp_minimize()`.

        bins : int, bins=20
            Number of bins to use for histograms on the diagonal.

        dimensions : list of str, default=None
            Labels of the dimension
            variables. `None` defaults to `space.dimensions[i].name`, or
            if also `None` to `['X_0', 'X_1', ..]`.

        plot_dims : list of str and int, default=None
            List of dimension names or dimension indices from the
            search-space dimensions to be included in the plot.
            If `None` then use all dimensions except constant ones
            from the search-space.

        size : float, default=2
            Height (in inches) of each facet.

        wspace : float, default=0.35
            The width of the padding between subplots, as a fraction of the
            average Axes width. Ignored if ``ax`` is provided.

        hspace : float, default=0.35
           The height of the padding between subplots, as a fraction of the
           average Axes height. Ignored if ``ax`` is provided.

        size : float, default=2
            Height (in inches) of each facet.

        cmap: str or Colormap, default = 'viridis'
            Color map for scatter plots. Passed directly to
            `plt.scatter()`

        ax: `Matplotlib.Axes`, default= None
            An axis object in which to plot the dependence plot.

        Returns
        -------
        ax : `Matplotlib.Axes`
            Matplotlib axis the plto was drawn in
        """
        space = result.space
        # Convert categoricals to integers, so we can ensure consistent ordering.
        # Assign indices to categories in the order they appear in the Dimension.
        # Matplotlib's categorical plotting functions are only present in v 2.1+,
        # and may order categoricals differently in different plots anyway.
        samples, minimum, iscat = _map_categories(space, result.x_iters, result.x)
        order = range(samples.shape[0])

        if plot_dims is None:
            # Get all dimensions.
            plot_dims = []
            for row in range(space.n_dims):
                if space.dimensions[row].is_constant:
                    continue
                plot_dims.append((row, space.dimensions[row]))
        else:
            plot_dims = space[plot_dims]
        # Number of search-space dimensions we are using.
        n_dims = len(plot_dims)
        if dimensions is not None:
            assert len(dimensions) == n_dims

        fig_kwargs = dict(figsize=(size * n_dims, size * n_dims))
        ax, axes = _make_subgrid(
            ax, n_dims, fig_kwargs_=fig_kwargs, wspace=wspace, hspace=hspace
        )

        for i in range(n_dims):
            for j in range(n_dims):
                if i == j:
                    index, dim = plot_dims[i]
                    if iscat[index]:
                        bins_ = len(dim.categories)
                    elif dim.prior == 'log-uniform':
                        low, high = space.bounds[index]
                        bins_ = np.logspace(np.log10(low), np.log10(high), bins)
                    else:
                        bins_ = bins
                    if n_dims == 1:
                        ax_ = axes
                    else:
                        ax_ = axes[i, i]
                    ax_.hist(
                        samples[:, index],
                        bins=bins_,
                        range=None if iscat[index] else dim.bounds,
                    )

                # lower triangle
                elif i > j:
                    index_i, dim_i = plot_dims[i]
                    index_j, dim_j = plot_dims[j]
                    ax_ = axes[i, j]
                    ax_.scatter(
                        samples[:, index_j],
                        samples[:, index_i],
                        c=order,
                        s=40,
                        lw=0.0,
                        cmap=cmap,
                    )
                    ax_.scatter(
                        minimum[index_j],
                        minimum[index_i],
                        c=['r'],
                        s=100,
                        lw=0.0,
                        marker='*',
                    )

        # Make various adjustments to the plots.
        _format_scatter_plot_axes(
            axes,
            space,
            ylabel="Number of samples",
            plot_dims=plot_dims,
            dim_labels=dimensions,
        )
        return ax

    def _get_ylim_diagonal(ax):
        """Get the min / max of the ylim for all diagonal plots. This is used in
        _adjust_fig() so the ylim is the same for all diagonal plots.

        Parameters
        ----------
        ax : `Matplotlib.Axes`
            2-dimensional matrix with Matplotlib Axes objects.

        Returns
        -------
        ylim_diagonal : tuple(int)
            The common min and max ylim for the diagonal plots.
        """

        # Number of search-space dimensions used in this plot.
        if isinstance(ax, (list, np.ndarray)):
            n_dims = len(ax)
            # Get ylim for all diagonal plots.
            ylim = [ax[row, row].get_ylim() for row in range(n_dims)]
        else:
            ylim = [ax.get_ylim()]

        # Separate into two lists with low and high ylim.
        ylim_lo, ylim_hi = zip(*ylim)

        # Min and max ylim for all diagonal plots.
        ylim_min = np.min(ylim_lo)
        ylim_max = np.max(ylim_hi)

        return ylim_min, ylim_max

    def partial_dependence_1D(space, model, i, samples, n_points=40):
        """Calculate the partial dependence for a single dimension.

        This uses the given model to calculate the average objective value
        for all the samples, where the given dimension is fixed at
        regular intervals between its bounds.

        This shows how the given dimension affects 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 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`.

        Returns
        -------
        xi : np.array
            The points at which the partial dependence was evaluated.

        yi : np.array
            The average value of the modelled objective function at
            each point `xi`.
        """
        # 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):
            """Helper-function to calculate the average predicted objective value for the
            given model, when setting the index'th dimension of the search-space to the
            value x, and then averaging over all samples."""
            rvs_ = np.array(samples)  # copy
            # We replace the values in the dimension that we want to keep
            # fixed
            rvs_[:, dim_locs[i] : dim_locs[i + 1]] = x
            # In case of `x_eval=None` rvs conists of random samples.
            # Calculating the mean of these samples is how partial dependence
            # is implemented.
            return np.mean(model.predict(rvs_))

        xi, xi_transformed = _evenly_sample(space.dimensions[i], n_points)
        # Calculate the partial dependence for all the points.
        yi = [_calc(x) for x in xi_transformed]

        return xi, yi

    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

    def plot_objective_2D(
        result,
        dimension_identifier1,
        dimension_identifier2,
        n_points=40,
        n_samples=250,
        levels=10,
        zscale='linear',
        sample_source='random',
        minimum='result',
        n_minimum_search=None,
        ax=None,
    ):
        """Create and return a Matplotlib figure and axes with a landscape contour-plot of
        the last fitted model of the search-space, overlaid with all the samples from the
        optimization results, for the two given dimensions of the search-space.

        This is similar to `plot_objective()` but only for 2 dimensions
        whose doc-string also has a more extensive explanation.

        Parameters
        ----------
        result : `OptimizeResult`
            The optimization results e.g. from calling `gp_minimize()`.

        dimension_identifier1 : str or int
            Name or index of a dimension in the search-space.

        dimension_identifier2 : str or int
            Name or index of a dimension in the search-space.

        n_samples : int, default=250
            Number of random samples used for estimating the contour-plot
            of the objective function.

        n_points : int, default=40
            Number of points along each dimension where the partial dependence
            is evaluated when generating the contour-plots.

        levels : int, default=10
            Number of levels to draw on the contour plot.

        zscale : str, default='linear'
            Scale to use for the z axis of the contour plots.
            Either 'log' or linear for all other choices.

        ax : `Matplotlib.Axes`, default: None
            When set, everything is plotted inside this axis.

        Returns
        -------
        ax : `Matplotlib.Axes`
            The Matplotlib Figure-object.
            For example, you can save the plot by calling
            `fig.savefig('file.png')`
        """

        # Get the search-space instance from the optimization results.
        space = result.space
        x_vals = _evaluate_min_params(result, minimum, n_minimum_search)
        if sample_source == "random":
            x_eval = None
            samples = space.transform(space.rvs(n_samples=n_samples))
        else:
            x_eval = _evaluate_min_params(result, sample_source, n_minimum_search)
            samples = space.transform([x_eval])
        x_samples, x_minimum, _ = _map_categories(space, result.x_iters, x_vals)
        # Get the dimension-object, its index in the search-space, and its name.
        index1, dimension1 = space[dimension_identifier1]
        index2, dimension2 = space[dimension_identifier2]

        # Get the samples from the optimization-log for the relevant dimensions.
        # samples1 = get_samples_dimension(result=result, index=index1)
        samples1 = x_samples[:, index1]
        samples2 = x_samples[:, index2]
        # samples2 = get_samples_dimension(result=result, index=index2)

        # Get the best-found samples for the relevant dimensions.
        best_sample1 = x_minimum[index1]
        best_sample2 = x_minimum[index2]

        # Get the last fitted model for the search-space.
        last_model = result.models[-1]

        # Estimate the objective function for these sampled points
        # using the last fitted model for the search-space.
        xi, yi, zi = partial_dependence_2D(
            space, last_model, index2, index1, samples, n_points=n_points
        )

        if ax is None:
            ax = plt.gca()

        # Scale for the z-axis of the contour-plot. Either Log or Linear (None).
        locator = LogLocator() if zscale == 'log' else None

        # Plot the contour-landscape for the objective function.
        ax.contourf(xi, yi, zi, levels, locator=locator, cmap='viridis_r')

        # Plot all the parameters that were sampled during optimization.
        # These are plotted as small black dots.
        ax.scatter(samples1, samples2, c='black', s=10, linewidths=1)

        # Plot the best parameters that were sampled during optimization.
        # These are plotted as a big red star.
        ax.scatter(best_sample1, best_sample2, c='red', s=50, linewidths=1, marker='*')

        # Use the dimension-names as the labels for the plot-axes.
        ax.set_xlabel(dimension1.name)
        ax.set_ylabel(dimension2.name)
        ax.autoscale(enable=True, axis='x', tight=True)
        ax.autoscale(enable=True, axis='y', tight=True)
        # Use log-scale on the x-axis?
        if dimension1.prior == 'log-uniform':
            ax.set_xscale('log')

        # Use log-scale on the y-axis?
        if dimension2.prior == 'log-uniform':
            ax.set_yscale('log')
        return ax

    def plot_histogram(result, dimension_identifier, bins=20, rotate_labels=0, ax=None):
        """Create and return a Matplotlib figure with a histogram of the samples from the
        optimization results, for a given dimension of the search-space.

        Parameters
        ----------
        result : `OptimizeResult`
            The optimization results e.g. from calling `gp_minimize()`.

        dimension_identifier : str or int
            Name or index of a dimension in the search-space.

        bins : int, bins=20
            Number of bins in the histogram.

        rotate_labels : int, rotate_labels=0
            Degree to rotate category-names on the x-axis.
            Only used for Categorical dimensions.

        Returns
        -------
        ax : `Matplotlib.Axes`
            The Matplotlib Axes-object.
        """

        # Get the search-space instance from the optimization results.
        space = result.space

        # Get the dimension-object.
        index, dimension = space[dimension_identifier]

        # Get the samples from the optimization-log for that particular dimension.
        samples = [x[index] for x in result.x_iters]

        if ax is None:
            ax = plt.gca()

        if isinstance(dimension, Categorical):
            # When the search-space dimension is Categorical, it means
            # that the possible values are strings. Matplotlib's histogram
            # does not support this, so we have to make a bar-plot instead.

            # NOTE: This only shows the categories that are in the samples.
            # So if a category was not sampled, it will not be shown here.

            # Count the number of occurrences of the string-categories.
            counter = Counter(samples)

            # The counter returns a dict where the keys are the category-names
            # and the values are the number of occurrences for each category.
            names = list(counter.keys())
            counts = list(counter.values())

            # Although Matplotlib's docs indicate that the bar() function
            # can take a list of strings for the x-axis, it doesn't appear to work.
            # So we hack it by creating a list of integers and setting the
            # tick-labels with the category-names instead.
            x = np.arange(len(counts))

            # Plot using bars.
            ax.bar(x, counts, tick_label=names)

            # Adjust the rotation of the category-names on the x-axis.
            ax.set_xticklabels(labels=names, rotation=rotate_labels)
        else:
            # Otherwise the search-space Dimension is either integer or float,
            # in which case the histogram can be plotted more easily.
            if dimension.prior == 'log-uniform':
                # Map the number of bins to a log-space for the dimension bounds.
                bins_mapped = np.logspace(*np.log10(dimension.bounds), bins)
            else:
                # Use the original number of bins.
                bins_mapped = bins
            # Plot the histogram.
            ax.hist(samples, bins=bins_mapped, range=dimension.bounds)

            # Use log-scale on the x-axis?
            if dimension.prior == 'log-uniform':
                ax.set_xscale('log')

        # Set the labels.
        ax.set_xlabel(dimension.name)
        ax.set_ylabel('Sample Count')

        return ax

    def _map_categories(space, points, minimum):
        """Map categorical values to integers in a set of points.

        Returns
        -------
        mapped_points : np.array, shape=points.shape
            A copy of `points` with categoricals replaced with their indices in
            the corresponding `Dimension`.

        mapped_minimum : np.array, shape (space.n_dims,)
            A copy of `minimum` with categoricals replaced with their indices in
            the corresponding `Dimension`.

        iscat : np.array, shape (space.n_dims,)
           Boolean array indicating whether dimension `i` in the `space` is
           categorical.
        """
        points = np.asarray(points, dtype=object)  # Allow slicing, preserve cats
        iscat = np.repeat(False, space.n_dims)
        min_ = np.zeros(space.n_dims)
        pts_ = np.zeros(points.shape)
        for i, dim in enumerate(space.dimensions):
            if isinstance(dim, Categorical):
                iscat[i] = True
                catmap = dict(zip(dim.categories, count()))
                pts_[:, i] = [catmap[cat] for cat in points[:, i]]
                min_[i] = catmap[minimum[i]]
            else:
                pts_[:, i] = points[:, i]
                min_[i] = minimum[i]
        return pts_, min_, iscat

    def _evenly_sample(dim, n_points):
        """Return `n_points` evenly spaced points from a Dimension.

        Parameters
        ----------
        dim : `Dimension`
            The Dimension to sample from.  Can be categorical; evenly-spaced
            category indices are chosen in order without replacement (result
            may be smaller than `n_points`).

        n_points : int
            The number of points to sample from `dim`.

        Returns
        -------
        xi : np.array
            The sampled points in the Dimension.  For Categorical
            dimensions, returns the index of the value in
            `dim.categories`.

        xi_transformed : np.array
            The transformed values of `xi`, for feeding to a model.
        """
        cats = np.array(getattr(dim, 'categories', []), dtype=object)
        if len(cats):  # Sample categoricals while maintaining order
            xi = np.linspace(0, len(cats) - 1, min(len(cats), n_points), dtype=int)
            xi_transformed = dim.transform(cats[xi])
        else:
            bounds = dim.bounds
            # XXX use linspace(*bounds, n_points) after python2 support ends
            xi = np.linspace(bounds[0], bounds[1], n_points)
            xi_transformed = dim.transform(xi)
        return xi, xi_transformed

    def _cat_format(dimension, x, _):
        """Categorical axis tick formatter function.

        Returns the name of category
        `x` in `dimension`.  Used with `matplotlib.ticker.FuncFormatter`.
        """
        if len(dimension.categories) > x:
            return str(dimension.categories[int(x)])
        else:
            return ''

    def _evaluate_min_params(
        result, params='result', n_minimum_search=None, random_state=None
    ):
        """Returns the minimum based on `params`"""
        x_vals = None
        if isinstance(params, str):
            if params == 'result':
                # Using the best observed result
                x_vals = result.x
            elif params == 'expected_minimum':
                if result.space.is_partly_categorical:
                    # space is also categorical
                    raise ValueError(
                        'expected_minimum does not support any' 'categorical values'
                    )
                # Do a gradient based minimum search using scipys own minimizer
                if n_minimum_search:
                    # If a value for
                    # expected_minimum_samples has been parsed
                    x_vals, _ = expected_minimum(
                        result, n_random_starts=n_minimum_search, random_state=random_state
                    )
                else:  # Use standard of 20 random starting points
                    x_vals, _ = expected_minimum(
                        result, n_random_starts=20, random_state=random_state
                    )
            elif params == 'expected_minimum_random':
                # Do a minimum search by evaluating the function with
                # n_samples sample values
                if n_minimum_search is not None:
                    # If a value for
                    # n_minimum_samples has been parsed
                    x_vals, _ = expected_minimum_random_sampling(
                        result, n_random_starts=n_minimum_search, random_state=random_state
                    )
                else:
                    # Use standard of 10^n_parameters. Note this
                    # becomes very slow for many parameters
                    n_minimum_search = 10 ** len(result.x)
                    x_vals, _ = expected_minimum_random_sampling(
                        result, n_random_starts=n_minimum_search, random_state=random_state
                    )
            else:
                raise ValueError(
                    'Argument ´eval_min_params´ must be a valid' 'string (´result´)'
                )
        elif isinstance(params, list):
            assert len(params) == len(result.x), (
                'Argument'
                '´eval_min_params´ of type list must have same length as'
                'number of features'
            )
            # Using defined x_values
            x_vals = params
        else:
            raise ValueError('Argument ´eval_min_params´ must' 'be a string or a list')
        return x_vals