diff --git a/CHANGELOG.md b/CHANGELOG.md index 92a52b2..30047a6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,19 @@ # Changelog +## [0.7.13] + +### Added +- Campaign.Explainer now added to PyTests +- Docstrings and typing to Explainer methods + +### Modified +- Fixed SHAP explainer analysis and visualization functions +- Changed SHAP visualization colors to use obsidian branding +- Moved sensitivity method from campaign.analysis to campaign.explainer + +### Removed +- Removed code chunks regarding unused optional inputs to PDP ICE function imported from SHAP GitHub + ## [0.7.12] ### Added - More informative docstrings for optimizer.bayesian, optimizer.predict, to explain choices of surrogate models and aq_funcs diff --git a/obsidian/campaign/analysis.py b/obsidian/campaign/analysis.py index c894d9a..07c663e 100644 --- a/obsidian/campaign/analysis.py +++ b/obsidian/campaign/analysis.py @@ -195,50 +195,3 @@ def calc_ofat_ranges(optimizer, threshold, X_ref, PI_range=0.7, cor = None return ofat_ranges, cor - - -def sensitivity(optimizer, - dx: float = 1e-6, - X_ref: pd.DataFrame | None = None) -> pd.DataFrame: - """ - Calculates the sensitivity of the surrogate model predictions with respect to each parameter in the X_space. - - Args: - optimizer (BayesianOptimizer): The optimizer object which contains a surrogate that has been fit to data - and can be used to make predictions. - dx (float, optional): The perturbation size for calculating the sensitivity. Defaults to ``1e-6``. - X_ref (pd.DataFrame | None, optional): The reference input values for calculating the sensitivity. - If None, the mean of X_space will be used as the reference. Defaults to ``None``. - - Returns: - pd.DataFrame: A DataFrame containing the sensitivity values for each parameter in X_space. - - Raises: - ValueError: If X_ref does not contain all parameters in optimizer.X_space or if X_ref is not a single row DataFrame. - """ - if X_ref is None: - X_ref = optimizer.X_space.mean() - else: - if not all(x in X_ref.columns for x in optimizer.X_space.X_names): - raise ValueError('X_ref must contain all parameters in X_space') - if X_ref.shape[0] != 1: - raise ValueError('X_ref must be a single row DataFrame') - - y_ref = optimizer.predict(X_ref) - - sens = {} - - # Only do positive perturbation, for simplicity - for param in optimizer.X_space: - base = param.unit_map(X_ref[param.name].values)[0] - # Space already mapped to (0,1), use absolute perturbation - dx_pos = np.array(base+dx).reshape(-1, 1) - X_sim = X_ref.copy() - X_sim[param.name] = param.unit_demap(dx_pos)[0] - y_sim = optimizer.predict(X_sim) - dydx = (y_sim - y_ref)/dx - sens[param.name] = dydx.to_dict('records')[0] - - df_sens = pd.DataFrame(sens).T[[y+' (pred)' for y in optimizer.y_names]] - - return df_sens diff --git a/obsidian/campaign/campaign.py b/obsidian/campaign/campaign.py index bc542a9..da7a9a3 100644 --- a/obsidian/campaign/campaign.py +++ b/obsidian/campaign/campaign.py @@ -29,6 +29,7 @@ class Campaign(): m_exp (int): The number of observations in campaign.data y (pd.Series): The response data in campaign.data f (pd.Series): The transformed response data + o (pd.Series): The objective function evaluated on f X (pd.DataFrame): The input features of campaign.data response_max (float | pd.Series): The maximum for each response target (Target | list[Target]): The target(s) for optimization. @@ -207,7 +208,7 @@ def f(self) -> pd.Series | pd.DataFrame: def o(self) -> pd.Series | pd.DataFrame: if self.objective: try: - x = self.X_space.encode(self.data[list(self.X_space.X_names)]).values + x = self.X_space.encode(self.X).values o = self.objective(torch.tensor(self.f.values).unsqueeze(0), X=torch.tensor(x)).squeeze(0) if o.ndim < 2: @@ -224,7 +225,7 @@ def X(self) -> pd.DataFrame: """ Feature columns of the training data """ - return self.data[self.X_space.X_names] + return self.data[list(self.X_space.X_names)] def save_state(self) -> dict: """ diff --git a/obsidian/campaign/explainer.py b/obsidian/campaign/explainer.py index c37ba0c..369954b 100644 --- a/obsidian/campaign/explainer.py +++ b/obsidian/campaign/explainer.py @@ -1,7 +1,5 @@ """Explainer Class: Surrogate Model Interpretation Methods""" -from .analysis import sensitivity - from obsidian.parameters import Param_Continuous, ParamSpace from obsidian.optimizer import Optimizer @@ -35,7 +33,7 @@ class Explainer(): """ def __init__(self, optimizer: Optimizer, - X_space: ParamSpace | None = None): + X_space: ParamSpace | None = None) -> None: if not optimizer.is_fit: raise ValueError('Surrogate model in optimizer is not fit to data. ') @@ -45,15 +43,15 @@ def __init__(self, self.shap = {} - def __repr__(self): + def __repr__(self) -> str: return f"Explainer(optimizer={self.optimizer})" @property - def optimizer(self): + def optimizer(self) -> Optimizer: """Explainer Optimizer object""" return self._optimizer - def set_optimizer(self, optimizer: Optimizer): + def set_optimizer(self, optimizer: Optimizer) -> None: """Sets the explainer optimizer""" self._optimizer = optimizer @@ -61,7 +59,7 @@ def shap_explain(self, responseid: int = 0, n: int = 100, X_ref: pd.DataFrame | None = None, - seed: int | None = None): + seed: int | None = None) -> None: """ Explain the parameter sensitivities using shap values. @@ -93,7 +91,7 @@ def shap_explain(self, y_name = self.optimizer.target[responseid].name - def f_preds(X): + def pred_func(X): # helper function for shap_explain if isinstance(X, np.ndarray): X = pd.DataFrame(X, columns=self.X_space.X_names) @@ -104,13 +102,14 @@ def f_preds(X): mu = y_pred[y_name+' (pred)'].values return mu - self.shap['f_preds'] = f_preds - self.shap['explainer'] = KernelExplainer(f_preds, X_ref) + self.shap['pred_func'] = pred_func + self.shap['explainer'] = KernelExplainer(pred_func, X_ref) self.shap['X_sample'] = self.X_space.unit_demap( pd.DataFrame(torch.rand(size=(n, self.X_space.n_dim)).numpy(), columns=X_ref.columns) ) - self.shap['values'] = self.shap['explainer'].shap_values(self.shap['X_sample'], seed=seed, l1_reg='aic') + self.shap['values'] = self.shap['explainer'].shap_values(self.shap['X_sample'], + seed=seed, l1_reg='aic') self.shap['explanation'] = Explanation(self.shap['values'], feature_names=X_ref.columns) return @@ -120,7 +119,11 @@ def shap_summary(self) -> Figure: if not self.shap: raise ValueError('shap explainer is not fit.') - fig = shap.summary_plot(self.shap['values'], self.shap['X_sample']) + fig = plt.figure() + shap.summary_plot(self.shap['values'], self.shap['X_sample'], + show=False) + plt.close(fig) + return fig def shap_summary_bar(self) -> Figure: @@ -129,64 +132,148 @@ def shap_summary_bar(self) -> Figure: raise ValueError('shap explainer is not fit.') fig = plt.figure() - shap.plots.bar(self.shap['explanation'], ax=fig.gca(), show=False) + shap.plots.bar(self.shap['explanation'], + ax=fig.gca(), show=False) plt.close(fig) return fig def shap_pdp_ice(self, - ind=0, # var name or index - ice_color_var=0, # var name or index + ind: int | tuple[int] = 0, + ice_color_var: int = 0, ace_opacity: float = 0.5, - ace_linewidth="auto" + npoints: int | None = None, ) -> Figure: - """SHAP Partial Dependence Plot with ICE""" + """ + SHAP Partial Dependence Plot with ICE + + Args: + ind (int): Index of the parameter to plot + ice_color_var (int): Index of the parameter to color the ICE lines + ace_opacity (float): Opacity of the ACE line + npoints (int, optional): Number of points for PDP x-axis. By default + will use ``100`` for 1D PDP and ``20`` for 2D PDP. + + Returns: + Matplotlib Figure of 1D or 2D PDP with ICE lines + + """ if not self.shap: raise ValueError('shap explainer is not fit.') - fig = partial_dependence( + fig, ax = partial_dependence( ind=ind, - model=self.f_preds, - data=self.X_sample, + model=self.shap['pred_func'], + data=self.shap['X_sample'], ice_color_var=ice_color_var, hist=False, ace_opacity=ace_opacity, - ace_linewidth=ace_linewidth, - show=False - ) + show=False, + npoints=npoints + ) + plt.close(fig) return fig def shap_single_point(self, - X_new, - X_ref=None): - """SHAP Pair-wise Marginal Explanations""" + X_new: pd.DataFrame | pd.Series, + X_ref=None) -> tuple[pd.DataFrame, Figure, Figure]: + """ + SHAP Pair-wise Marginal Explanations + + Args: + X_new (pd.DataFrame | pd.Series): New data point to explain + X_ref (pd.DataFrame | pd.Series, optional): Reference data point + for shap values. Default uses ``X_space.mean()`` + + Returns: + pd.DataFrame: DataFrame containing SHAP values for the new data point + Figure: Matplotlib Figure for SHAP values + Figure: Matplotlib Figure for SHAP summary plot + + """ if not self.shap: raise ValueError('shap explainer is not fit.') + + if isinstance(X_new, pd.Series): + X_new = X_new.copy().to_frame().T + if isinstance(X_ref, pd.Series): + X_ref = X_ref.copy().to_frame().T + + if not list(X_new.columns) == list(self.optimizer.X_space.X_names): + raise ValueError('X_new must contain all parameters in X_space') if X_ref is None: - if self.shap_explainer is None: - raise ValueError('shap explainer is not fit. ') - return - shap_value_new = self.shap_explainer.shap_values(X_new) - expected_value = self.shap_explainer.expected_value + shap_value_new = self.shap['explainer'].shap_values(X_new).squeeze() + expected_value = self.shap['explainer'].expected_value else: + if not list(X_ref.columns) == list(self.optimizer.X_space.X_names): + raise ValueError('X_ref must contain all parameters in X_space') + # if another reference point is input, need to re-fit another explainer - explainer = shap.KernelExplainer(self.f_preds, X_ref) - shap_value_new = explainer.shap_values(X_new) + explainer = shap.KernelExplainer(self.shap['pred_func'], X_ref) + shap_value_new = explainer.shap_values(X_new).squeeze() expected_value = explainer.expected_value - shap_value_new = np.squeeze(shap_value_new) df_shap_value_new = pd.DataFrame([shap_value_new], columns=self.X_space.X_names) - fig1, fig2 = one_shap_value(shap_value_new, expected_value, self.X_space.X_names) + fig_bar, fig_line = one_shap_value(shap_value_new, expected_value, self.X_space.X_names) - return df_shap_value_new, fig1, fig2 + return df_shap_value_new, fig_bar, fig_line + + def sensitivity(self, + dx: float = 1e-6, + X_ref: pd.DataFrame | pd.Series | None = None) -> pd.DataFrame: + """ + Calculates the local sensitivity of the surrogate model predictions with + respect to each parameter in the X_space. + + Args: + optimizer (BayesianOptimizer): The optimizer object which contains a surrogate + that has been fit to data + and can be used to make predictions. + dx (float, optional): The perturbation size for calculating the sensitivity. + Defaults to ``1e-6``. + X_ref (pd.DataFrame | pd.Series | None, optional): The reference input values for + calculating the sensitivity. If None, the mean of X_space will be used as the + reference. Defaults to ``None``. - def cal_sensitivity(self, - dx: float = 1e-6, - X_ref: pd.DataFrame | None = None) -> pd.DataFrame: - """Local parameter sensitivity analysis""" - df_sens = sensitivity(self.optimizer, dx=dx, X_ref=X_ref) + Returns: + pd.DataFrame: A DataFrame containing the sensitivity values for each parameter + in X_space. + + Raises: + ValueError: If X_ref does not contain all parameters in optimizer.X_space or if + X_ref is not a single row DataFrame. + + """ + + if isinstance(X_ref, pd.Series): + X_ref = X_ref.copy().to_frame().T + + if X_ref is None: + X_ref = self.optimizer.X_space.mean() + else: + if not all(x in X_ref.columns for x in self.optimizer.X_space.X_names): + raise ValueError('X_ref must contain all parameters in X_space') + if X_ref.shape[0] != 1: + raise ValueError('X_ref must be a single row DataFrame') + + y_ref = self.optimizer.predict(X_ref) + + sens = {} + + # Only do positive perturbation, for simplicity + for param in self.optimizer.X_space: + base = param.unit_map(X_ref[param.name].values)[0] + # Space already mapped to (0,1), use absolute perturbation + dx_pos = np.array(base+dx).reshape(-1, 1) + X_sim = X_ref.copy() + X_sim[param.name] = param.unit_demap(dx_pos)[0] + y_sim = self.optimizer.predict(X_sim) + dydx = (y_sim - y_ref)/dx + sens[param.name] = dydx.to_dict('records')[0] + + df_sens = pd.DataFrame(sens).T[[y+' (pred)' for y in self.optimizer.y_names]] return df_sens diff --git a/obsidian/plotting/shap.py b/obsidian/plotting/shap.py index d5df883..6f11785 100644 --- a/obsidian/plotting/shap.py +++ b/obsidian/plotting/shap.py @@ -1,105 +1,140 @@ +"""Custom plots for SHAP analysis visualization""" + +from .branding import obsidian_cm, obsidian_colors + +from shap.plots._partial_dependence import compute_bounds +from shap.utils import convert_name + import matplotlib.pyplot as plt +from matplotlib.axes import Axes +from matplotlib.figure import Figure + import numpy as np import pandas as pd -from shap import Explanation -from shap.plots.colors import blue_rgb, light_blue_rgb, red_blue_transparent, red_rgb -from shap.utils import convert_name -from shap.plots._partial_dependence import compute_bounds - -from matplotlib.pyplot import get_cmap +from typing import Callable -def one_shap_value(shap_value_new, expected_value, X_names): +def one_shap_value(shap_value_new: np.ndarray, + expected_value: float, + X_names: list[str]) -> tuple[Figure, Figure]: """ Visualize the shap values of one data point + + Args: + shap_value_new (np.ndarray): The SHAP values of a single data point + to be compared to a reference point. + expected_value (float): The expected value at the reference point. + X_names (list[str]): The names of the features. + + Returns: + Figure: The bar plot of SHAP values for the single data point. + Figure: The line plot of cumulative SHAP values for the data point in + comparison to the reference point. + """ - cumulative_shap = np.cumsum(shap_value_new) + # First figure = Bar plot, just SHAP values of the new point pred_new = expected_value + np.sum(shap_value_new) - fig1 = plt.figure(figsize=(8, 6)) - plt.bar(range(len(shap_value_new)), shap_value_new) + fig_bar = plt.figure(figsize=(8, 6)) + plt.bar(range(len(shap_value_new)), shap_value_new, color=obsidian_colors.primary.teal) plt.xticks(range(len(shap_value_new)), X_names) plt.ylabel('SHAP Value') plt.xlabel('Feature') plt.title('SHAP Values for Single Datapoint') plt.axhline(0, color='grey', linestyle='--') - plt.close(fig1) + plt.close(fig_bar) - fig2 = plt.figure(figsize=(8, 6)) + # Second figure = Line plot, cumulative SHAP values from new point to reference + fig_line = plt.figure(figsize=(8, 6)) + cumulative_shap = np.cumsum(shap_value_new) y = cumulative_shap+expected_value n = len(cumulative_shap) - plt.plot(range(n), y, color = 'steelblue') - plt.scatter(range(n), y, color='blue') + plt.plot(range(n), y, color=obsidian_colors.primary.teal) + plt.scatter(range(n), y, color=obsidian_colors.secondary.light_teal) plt.vlines(0, expected_value, y[0], colors='steelblue', linestyles='solid') - plt.scatter(0, expected_value, color='purple') - plt.scatter(n-1, y[-1], color='pink') + plt.scatter(0, expected_value, color=obsidian_colors.accent.vista_blue) + plt.scatter(n-1, y[-1], color=obsidian_colors.secondary.blue) plt.xticks(range(n), X_names) plt.xlabel('Feature') plt.ylabel('Cumulative SHAP Value') plt.title('Cumulative SHAP Values for Single Datapoint') - plt.axhline(expected_value, color='purple', linestyle='--',label='Reference') - plt.axhline(pred_new, color='pink', linestyle='--',label='New Data') - plt.ylim([min(min(y),expected_value)*0.98, max(max(y),expected_value)*1.02]) + plt.axhline(expected_value, color=obsidian_colors.accent.vista_blue, + linestyle='--', label='Reference') + plt.axhline(pred_new, color=obsidian_colors.secondary.blue, + linestyle='--', label='New Data') + plt.ylim([min(min(y), expected_value)*0.98, max(max(y), expected_value)*1.02]) plt.legend() - plt.close(fig2) + plt.close(fig_line) - return fig1, fig2 - - -def partial_dependence(ind, model, data, - ice_color_var=None, - xmin="percentile(0)", xmax="percentile(100)", - npoints=None, feature_names=None, hist=True, model_expected_value=False, - feature_expected_value=False, shap_values=None, - ylabel=None, ice=True, ace_opacity=1, pd_opacity=1, pd_linewidth=2, - ace_linewidth='auto', ax=None, show=True): + return fig_bar, fig_line + + +def partial_dependence(ind: int | tuple[int], + model: Callable, + data: pd.DataFrame, + ice_color_var: int | str | None = 0, + xmin: str | tuple[float] | float = "percentile(0)", + xmax: str | tuple[float] | float = "percentile(100)", + npoints: int | None = None, + feature_names: list[str] | None = None, + hist: bool = False, + ylabel: str | None = None, + ice: bool = True, + ace_opacity: float = 1, + pd_opacity: float = 1, + pd_linewidth: float = 2, + ace_linewidth: str | float = 'auto', + ax: Axes | None = None, + show: bool = True) -> Figure: """ Calculates and plots the partial dependence of a feature or a pair of features on the model's output. - This function is revised from the partial_dependence_plot function in shap package, - in order to color the ICE curves by certain feature for checking interaction between features. + This function is revised from the partial_dependence_plot function in shap package, + in order to color the ICE curves by certain feature for checking interaction between features. Ref: https://github.com/shap/shap/blob/master/shap/plots/_partial_dependence.py Args: - ind (int or tuple): The index or indices of the feature(s) to calculate the partial dependence for. - model: The model used for prediction. - data: The input data used for prediction. - ice_color_var: The index of the feature used for coloring the ICE lines (for 1D partial dependence plot). - xmin (str or tuple or float): The minimum value(s) for the feature(s) range. - xmax (str or tuple or float): The maximum value(s) for the feature(s) range. - npoints (int): The number of points to sample within the feature(s) range. - feature_names (list): The names of the features. - hist (bool): Whether to plot the histogram of the feature(s). - model_expected_value (bool or float): Whether to plot the model's expected value line. - feature_expected_value (bool): Whether to plot the feature's expected value line. - shap_values: The SHAP values used for plotting. - ylabel (str): The label for the y-axis. - ice (bool): Whether to plot the Individual Conditional Expectation (ICE) lines. - ace_opacity (float): The opacity of the ACE lines. - pd_opacity (float): The opacity of the PDP line. - pd_linewidth (float): The linewidth of the PDP line. - ace_linewidth (float or str): The linewidth of the ACE lines. Set to 'auto' for automatic calculation. - ax: The matplotlib axis to plot on. - show (bool): Whether to show the plot. + ind (int | tuple): The index or indices of the feature(s) to calculate + the partial dependence for. + model (Callable): The model used for prediction. + data (pd.DataFrame): The input data used for prediction. + ice_color_var (int | str, optional): The index of the feature used for coloring + the ICE lines (for 1D partial dependence plot). Default is ``0``. + xmin (str | tuple | float, optional): The minimum value(s) for the feature(s) range. + Default is ``"percentile(0)"``. + xmax (str | tuple | float): The maximum value(s) for the feature(s) range. + Default is ``"percentile(100)"``. + npoints (int, optional): The number of points to sample within the feature(s) range. + By default, will use ``100`` points for 1D PDP and ``20`` points for 2D PDP. + feature_names (list[str], optional): The names of the features. Will default to + the names of the DataFrame columns. + hist (bool, optional): Whether to plot the histogram of the feature(s). Default + is ``False``. + ylabel (str, optional): The label for the y-axis. Default is ``None``. + ice (bool, optional): Whether to plot the Individual Conditional Expectation (ICE) lines. + Default is ``True``. + ace_opacity (float, optional): The opacity of the ACE lines. Default is ``1``. + pd_opacity (float, optional): The opacity of the PDP line. Default is ``1``. + pd_linewidth (float, optional): The linewidth of the PDP line. Default is ``2``. + ace_linewidth (float | str, optional): The linewidth of the ACE lines. Default is ``'auto'`` + for automatic calculation. + ax (Axes, optional): The matplotlib axis to plot on. By default will attach to Figure.gca(). + show (bool, optional): Whether to show the plot. Default is ``True``. Returns: tuple: A tuple containing the matplotlib figure and axis objects if `show` is False, otherwise None. """ - if isinstance(data, Explanation): - features = data.data - shap_values = data - else: - features = data + features = data - # convert from DataFrames if we got any + # Convert from DataFrame if used use_dataframe = False if isinstance(features, pd.DataFrame): if feature_names is None: @@ -110,7 +145,7 @@ def partial_dependence(ind, model, data, if feature_names is None: feature_names = ["Feature %d" % i for i in range(features.shape[1])] - # this is for a 1D partial dependence plot + # 1D PDP if type(ind) is not tuple: ind = convert_name(ind, None, feature_names) xv = features[:, ind] @@ -127,10 +162,6 @@ def partial_dependence(ind, model, data, ice_vals[i, :] = model(pd.DataFrame(features_tmp, columns=feature_names)) else: ice_vals[i, :] = model(features_tmp) - # if linewidth is None: - # linewidth = 1 - # if opacity is None: - # opacity = 0.5 features_tmp = features.copy() vals = np.zeros(npoints) @@ -148,7 +179,6 @@ def partial_dependence(ind, model, data, fig = plt.gcf() ax1 = plt.gca() - # fig, ax1 = plt.subplots(figsize) ax2 = ax1.twinx() # the histogram of the data @@ -159,32 +189,33 @@ def partial_dependence(ind, model, data, # ice line plot if ice: if ace_linewidth == "auto": - ace_linewidth = min(1, 50/ice_vals.shape[1]) # pylint: disable=unsubscriptable-object - # -------------------------- Revised + ace_linewidth = min(1, 50/ice_vals.shape[1]) + if ice_color_var is None: - ax1.plot(xs, ice_vals, color=light_blue_rgb, linewidth=ace_linewidth, alpha=ace_opacity) + ax1.plot(xs, ice_vals, color=obsidian_colors.secondary.light_teal, + linewidth=ace_linewidth, alpha=ace_opacity) else: if not isinstance(ice_color_var, int): ice_color_var = convert_name(ice_color_var, None, feature_names) - colormap = get_cmap('coolwarm') + colormap = obsidian_cm.obsidian_viridis color_vals = features[:, ice_color_var] colorbar_min = color_vals.min() colorbar_max = color_vals.max() color_vals = (color_vals - colorbar_min)/(colorbar_max-colorbar_min) - for i in range(ice_vals.shape[0]): + for i in range(ice_vals.shape[1]): ax1.plot(xs, ice_vals[:, i], color=colormap(color_vals[i]), linewidth=ace_linewidth, alpha=ace_opacity) cbar = plt.colorbar(plt.cm.ScalarMappable(cmap=colormap, - norm=plt.Normalize( - vmin=colorbar_min, vmax=colorbar_max)), ax=ax1) + norm=plt.Normalize( + vmin=colorbar_min, vmax=colorbar_max)), + ax=ax1) cbar.set_label('Color by ' + feature_names[ice_color_var]) - # -------------------------- # the line plot ax1.plot(xs, vals, color="black", linewidth=pd_linewidth, alpha=pd_opacity) # Revised: PDP line from blue_rgb to black - ax2.set_ylim(0, features.shape[0]) # ax2.get_ylim()[0], ax2.get_ylim()[1] * 4) + ax2.set_ylim(0, features.shape[0]) ax1.set_xlabel(feature_names[ind], fontsize=13) if ylabel is None: if not ice: @@ -206,66 +237,12 @@ def partial_dependence(ind, model, data, ax2.spines['left'].set_visible(False) ax2.spines['bottom'].set_visible(False) - if feature_expected_value is not False: - ax3 = ax2.twiny() - ax3.set_xlim(xmin, xmax) - mval = xv.mean() - ax3.set_xticks([mval]) - ax3.set_xticklabels(["E["+str(feature_names[ind])+"]"]) - ax3.spines['right'].set_visible(False) - ax3.spines['top'].set_visible(False) - ax3.tick_params(length=0, labelsize=11) - ax1.axvline(mval, color="#999999", zorder=-1, linestyle="--", linewidth=1) - - if model_expected_value is not False or shap_values is not None: - if model_expected_value is True: - if use_dataframe: - model_expected_value = model(pd.DataFrame(features, columns=feature_names)).mean() - else: - model_expected_value = model(features).mean() - else: - model_expected_value = shap_values.base_values - ymin, ymax = ax1.get_ylim() - ax4 = ax2.twinx() - ax4.set_ylim(ymin, ymax) - ax4.set_yticks([model_expected_value]) - ax4.set_yticklabels(["E[f(x)]"]) - ax4.spines['right'].set_visible(False) - ax4.spines['top'].set_visible(False) - ax4.tick_params(length=0, labelsize=11) - ax1.axhline(model_expected_value, color="#999999", zorder=-1, linestyle="--", linewidth=1) - - if shap_values is not None: - # vals = shap_values.values[:, ind] - # if shap_value_features is None: - # shap_value_features = features - # assert shap_values.shape == features.shape - # #sample_ind = 18 - # vals = shap_values[:, ind] - # if type(model_expected_value) is bool: - # if use_dataframe: - # model_expected_value = model(pd.DataFrame(features, columns=feature_names)).mean() - # else: - # model_expected_value = model(features).mean() - # if isinstance(shap_value_features, pd.DataFrame): - # shap_value_features = shap_value_features.values - markerline, stemlines, _ = ax1.stem( - shap_values.data[:, ind], shap_values.base_values + shap_values.values[:, ind], - bottom=shap_values.base_values, - markerfmt="o", basefmt=" ", use_line_collection=True - ) - stemlines.set_edgecolors([red_rgb if v > 0 else blue_rgb for v in vals]) - plt.setp(stemlines, 'zorder', -1) - plt.setp(stemlines, 'linewidth', 2) - plt.setp(markerline, 'color', "black") - plt.setp(markerline, 'markersize', 4) - if show: plt.show() else: return fig, ax1 - # this is for a 2D partial dependence plot + # 2D PDP else: ind0 = convert_name(ind[0], None, feature_names) ind1 = convert_name(ind[1], None, feature_names) @@ -298,12 +275,7 @@ def partial_dependence(ind, model, data, fig = plt.figure() ax = fig.add_subplot(111, projection='3d') -# x = y = np.arange(-3.0, 3.0, 0.05) -# X, Y = np.meshgrid(x, y) -# zs = np.array(fun(np.ravel(X), np.ravel(Y))) -# Z = zs.reshape(X.shape) - - ax.plot_surface(x0, x1, vals, cmap=red_blue_transparent) + ax.plot_surface(x0, x1, vals, cmap=obsidian_cm.obsidian_viridis) ax.set_xlabel(feature_names[ind0], fontsize=13) ax.set_ylabel(feature_names[ind1], fontsize=13) diff --git a/obsidian/tests/test_campaign.py b/obsidian/tests/test_campaign.py index d5dac53..37efa79 100644 --- a/obsidian/tests/test_campaign.py +++ b/obsidian/tests/test_campaign.py @@ -3,7 +3,7 @@ from obsidian.parameters import Target from obsidian.experiment import Simulator from obsidian.experiment.benchmark import two_leaves, shifted_parab -from obsidian.campaign import Campaign +from obsidian.campaign import Campaign, Explainer from obsidian.objectives import Identity_Objective, Scalar_WeightedNorm, Feature_Objective, \ Objective_Sequence, Utopian_Distance, Index_Objective, Bounded_Target @@ -72,6 +72,25 @@ def test_campaign_objectives(obj): campaign2.save_state() campaign2.__repr__() - + +def test_explain(): + exp = Explainer(campaign.optimizer) + exp.shap_explain(n=50) + + exp.shap_summary() + fig = exp.shap_summary_bar() + exp.shap_pdp_ice(ind=0, ice_color_var=None, npoints=10) + exp.shap_pdp_ice(ind=0, npoints=10) + exp.shap_pdp_ice(ind=(0, 1), npoints=5) + + X_new = campaign.X.iloc[0, :] + X_ref = campaign.X.loc[1, :] + df_shap_value_new, fig_bar, fig_line = exp.shap_single_point(X_new) + df_shap_value_new, fig_bar, fig_line = exp.shap_single_point(X_new, X_ref=X_ref) + + df_sens = exp.sensitivity() + df_sens = exp.sensitivity(X_ref=X_ref) + + if __name__ == '__main__': pytest.main([__file__, '-m', 'not slow'])