diff --git a/climada/test/test_util_calibrate.py b/climada/test/test_util_calibrate.py index 07a06c7737..61171e1aa7 100644 --- a/climada/test/test_util_calibrate.py +++ b/climada/test/test_util_calibrate.py @@ -11,7 +11,13 @@ from climada.entity import ImpactFuncSet, ImpactFunc -from climada.util.calibrate import Input, ScipyMinimizeOptimizer, BayesianOptimizer, OutputEvaluator +from climada.util.calibrate import ( + Input, + ScipyMinimizeOptimizer, + BayesianOptimizer, + OutputEvaluator, + BayesianOptimizerOutputEvaluator, +) from climada.util.calibrate.test.test_base import hazard, exposure @@ -213,6 +219,8 @@ def test_plots(self): output_eval.impf_set.plot() output_eval.plot_at_event() output_eval.plot_at_region() + output_eval.plot_event_region_heatmap() + + output_eval = BayesianOptimizerOutputEvaluator(self.input, output) ax = output_eval.plot_impf_variability() self.assertIsInstance(ax, Axes) - output_eval.plot_event_region_heatmap() diff --git a/climada/util/calibrate/__init__.py b/climada/util/calibrate/__init__.py index d1aa95d7dc..0c7100e023 100644 --- a/climada/util/calibrate/__init__.py +++ b/climada/util/calibrate/__init__.py @@ -1,6 +1,5 @@ """Impact function calibration module""" from .base import Input, OutputEvaluator -from .bayesian_optimizer import BayesianOptimizer +from .bayesian_optimizer import BayesianOptimizer, BayesianOptimizerOutputEvaluator from .scipy_optimizer import ScipyMinimizeOptimizer -from .func import rmse, rmsf, impact_at_reg diff --git a/climada/util/calibrate/base.py b/climada/util/calibrate/base.py index a2650a775f..540c75e72f 100644 --- a/climada/util/calibrate/base.py +++ b/climada/util/calibrate/base.py @@ -201,142 +201,6 @@ def __post_init__(self): ).impact(assign_centroids=True, save_mat=True) self._impact_label = f"Impact [{self.input.exposure.value_unit}]" - def plot_impf_variability( - self, - cost_func_diff: float = 0.1, - p_space_df: Optional[pd.DataFrame] = None, - plot_haz: bool = True, - plot_impf_kws: Optional[dict] = None, - plot_hist_kws: Optional[dict] = None, - ): - """Plot impact function variability with parameter combinations of - almost equal cost function values - - Args: - cost_func_diff (float, optional): Max deviation from optimal cost - function value (as fraction). Defaults to 0.1 (i.e. 10%). - p_space_df (pd.DataFrame, optional): parameter space. Defaults to None. - plot_haz (bool, optional): Whether or not to plot hazard intensity - distibution. Defaults to False. - plot_impf_kws (dict, optional): Keyword arguments for impact - function plot. Defaults to None. - plot_hist_kws (dict, optional): Keyword arguments for hazard - intensity distribution plot. Defaults to None. - """ - - # Initialize plot keyword arguments - if plot_impf_kws is None: - plot_impf_kws = {} - if plot_hist_kws is None: - plot_hist_kws = {} - - # Retrieve hazard type and parameter space - haz_type = self.input.hazard.haz_type - if p_space_df is None: - # Assert that self.output has the p_space_to_dataframe() method, - # which is defined for the BayesianOptimizerOutput class - if not hasattr(self.output, "p_space_to_dataframe"): - raise TypeError( - "To derive the full impact function parameter space, " - "plot_impf_variability() requires BayesianOptimizerOutput " - "as OutputEvaluator.output attribute, which provides the " - "method p_space_to_dataframe()." - ) - p_space_df = self.output.p_space_to_dataframe() - - # Retrieve parameters of impact functions with cost function values - # within 'cost_func_diff' % of the best estimate - params_within_range = p_space_df["Parameters"] - plot_space_label = "Parameter space" - if cost_func_diff is not None: - max_cost_func_val = p_space_df["Calibration", "Cost Function"].min() * ( - 1 + cost_func_diff - ) - params_within_range = params_within_range.loc[ - p_space_df["Calibration", "Cost Function"] <= max_cost_func_val - ] - plot_space_label = ( - f"within {int(cost_func_diff*100)} percent " f"of best fit" - ) - - # Set plot defaults - color = plot_impf_kws.pop("color", "tab:blue") - lw = plot_impf_kws.pop("lw", 2) - zorder = plot_impf_kws.pop("zorder", 3) - label = plot_impf_kws.pop("label", "best fit") - - # get number of impact functions and create a plot for each - n_impf = len(self.impf_set.get_func(haz_type=haz_type)) - axes = [] - - for impf_idx in range(n_impf): - _, ax = plt.subplots() - - # Plot best-fit impact function - best_impf = self.impf_set.get_func(haz_type=haz_type)[impf_idx] - ax.plot( - best_impf.intensity, - best_impf.mdd * best_impf.paa * 100, - color=color, - lw=lw, - zorder=zorder, - label=label, - **plot_impf_kws, - ) - - # Plot all impact functions within 'cost_func_diff' % of best estimate - for row in range(params_within_range.shape[0]): - label_temp = plot_space_label if row == 0 else None - - sel_params = params_within_range.iloc[row, :].to_dict() - temp_impf_set = self.input.impact_func_creator(**sel_params) - temp_impf = temp_impf_set.get_func(haz_type=haz_type)[impf_idx] - - ax.plot( - temp_impf.intensity, - temp_impf.mdd * temp_impf.paa * 100, - color="grey", - alpha=0.4, - label=label_temp, - ) - - # Plot hazard intensity value distributions - if plot_haz: - haz_vals = self.input.hazard.intensity[ - :, self.input.exposure.gdf[f"centr_{haz_type}"] - ] - - # Plot defaults - color_hist = plot_hist_kws.pop("color", "tab:orange") - alpha_hist = plot_hist_kws.pop("alpha", 0.3) - - ax2 = ax.twinx() - ax2.hist( - haz_vals.data, - bins=40, - color=color_hist, - alpha=alpha_hist, - label="Hazard intensity\noccurence", - ) - ax2.set(ylabel="Hazard intensity occurence (#Exposure points)") - ax.axvline( - x=haz_vals.max(), label="Maximum hazard value", color="tab:orange" - ) - ax2.legend(loc="lower right") - - ax.set( - xlabel=f"Intensity ({self.input.hazard.units})", - ylabel="Mean Damage Ratio (MDR) in %", - xlim=(min(best_impf.intensity), max(best_impf.intensity)), - ) - ax.legend() - axes.append(ax) - - if n_impf > 1: - return axes - - return ax - def plot_at_event( self, data_transf: Callable[[pd.DataFrame], pd.DataFrame] = lambda x: x, diff --git a/climada/util/calibrate/bayesian_optimizer.py b/climada/util/calibrate/bayesian_optimizer.py index 057ac4e467..913ddc9bac 100644 --- a/climada/util/calibrate/bayesian_optimizer.py +++ b/climada/util/calibrate/bayesian_optimizer.py @@ -6,129 +6,12 @@ from itertools import combinations, repeat import pandas as pd +import matplotlib.pyplot as plt import matplotlib.axes as maxes from bayes_opt import BayesianOptimization from bayes_opt.target_space import TargetSpace -from .base import Output, Optimizer - - -@dataclass -class BayesianOptimizer(Optimizer): - """An optimization using ``bayes_opt.BayesianOptimization`` - - This optimizer reports the target function value for each parameter set and - *maximizes* that value. Therefore, a higher target function value is better. - The cost function, however, is still minimized: The target function is defined as - the inverse of the cost function. - - For details on the underlying optimizer, see - https://github.com/bayesian-optimization/BayesianOptimization. - - Parameters - ---------- - input : Input - The input data for this optimizer. See the Notes below for input requirements. - verbose : int, optional - Verbosity of the optimizer output. Defaults to 1. - random_state : int, optional - Seed for initializing the random number generator. Defaults to 1. - allow_duplicate_points : bool, optional - Allow the optimizer to sample the same points in parameter space multiple times. - This may happen if the parameter space is tightly bound or constrained. Defaults - to ``True``. - bayes_opt_kwds : dict - Additional keyword arguments passed to the ``BayesianOptimization`` constructor. - - Notes - ----- - The following requirements apply to the parameters of :py:class:`Input` when using - this class: - - bounds - Setting ``bounds`` in the ``Input`` is required because the optimizer first - "explores" the bound parameter space and then narrows its search to regions - where the cost function is low. - constraints - Must be an instance of ``scipy.minimize.LinearConstraint`` or - ``scipy.minimize.NonlinearConstraint``. See - https://github.com/bayesian-optimization/BayesianOptimization/blob/master/examples/constraints.ipynb - for further information. Supplying contraints is optional. - - Attributes - ---------- - optimizer : bayes_opt.BayesianOptimization - The optimizer instance of this class. - """ - - verbose: InitVar[int] = 1 - random_state: InitVar[int] = 1 - allow_duplicate_points: InitVar[bool] = True - bayes_opt_kwds: InitVar[Optional[Mapping[str, Any]]] = None - - def __post_init__( - self, verbose, random_state, allow_duplicate_points, bayes_opt_kwds - ): - """Create optimizer""" - if bayes_opt_kwds is None: - bayes_opt_kwds = {} - - if self.input.bounds is None: - raise ValueError("Input.bounds is required for this optimizer") - - self.optimizer = BayesianOptimization( - f=self._opt_func, - pbounds=self.input.bounds, - constraint=self.input.constraints, - verbose=verbose, - random_state=random_state, - allow_duplicate_points=allow_duplicate_points, - **bayes_opt_kwds, - ) - - def _target_func(self, true: pd.DataFrame, predicted: pd.DataFrame) -> Number: - """Invert the cost function because BayesianOptimization maximizes the target""" - return -self.input.cost_func(true, predicted) - - def run(self, **opt_kwargs): - """Execute the optimization - - ``BayesianOptimization`` *maximizes* a target function. Therefore, this class - inverts the cost function and used that as target function. The cost function is - still minimized. - - Parameters - ---------- - init_points : int, optional - Number of initial samples taken from the parameter space. Defaults to 10^N, - where N is the number of parameters. - n_iter : int, optional - Number of iteration steps after initial sampling. Defaults to 10^N, where N - is the number of parameters. - opt_kwargs - Further keyword arguments passed to ``BayesianOptimization.maximize``. - - Returns - ------- - output : BayesianOptimizerOutput - Optimization output. :py:attr:`BayesianOptimizerOutput.p_space` stores data - on the sampled parameter space. - """ - # Retrieve parameters - num_params = len(self.input.bounds) - init_points = opt_kwargs.pop("init_points", 10**num_params) - n_iter = opt_kwargs.pop("n_iter", 10**num_params) - - # Run optimizer - self.optimizer.maximize(init_points=init_points, n_iter=n_iter, **opt_kwargs) - - # Return output - opt = self.optimizer.max - return BayesianOptimizerOutput( - params=opt["params"], - target=opt["target"], - p_space=self.optimizer.space, - ) +from .base import Output, Optimizer, OutputEvaluator @dataclass @@ -184,7 +67,7 @@ def plot_p_space( min_def: Optional[Union[str, Tuple[str, str]]] = "Cost Function", min_fmt: str = "x", min_color: str = "r", - **plot_kwargs + **plot_kwargs, ) -> Union[maxes.Axes, List[maxes.Axes]]: """Plot the parameter space as scatter plot(s) @@ -273,3 +156,275 @@ def plot_single(x, y): # Iterate over parameter combinations return [plot_single(p_first, p_second) for p_first, p_second in iterable] + + +@dataclass +class BayesianOptimizer(Optimizer): + """An optimization using ``bayes_opt.BayesianOptimization`` + + This optimizer reports the target function value for each parameter set and + *maximizes* that value. Therefore, a higher target function value is better. + The cost function, however, is still minimized: The target function is defined as + the inverse of the cost function. + + For details on the underlying optimizer, see + https://github.com/bayesian-optimization/BayesianOptimization. + + Parameters + ---------- + input : Input + The input data for this optimizer. See the Notes below for input requirements. + verbose : int, optional + Verbosity of the optimizer output. Defaults to 1. + random_state : int, optional + Seed for initializing the random number generator. Defaults to 1. + allow_duplicate_points : bool, optional + Allow the optimizer to sample the same points in parameter space multiple times. + This may happen if the parameter space is tightly bound or constrained. Defaults + to ``True``. + bayes_opt_kwds : dict + Additional keyword arguments passed to the ``BayesianOptimization`` constructor. + + Notes + ----- + The following requirements apply to the parameters of :py:class:`Input` when using + this class: + + bounds + Setting ``bounds`` in the ``Input`` is required because the optimizer first + "explores" the bound parameter space and then narrows its search to regions + where the cost function is low. + constraints + Must be an instance of ``scipy.minimize.LinearConstraint`` or + ``scipy.minimize.NonlinearConstraint``. See + https://github.com/bayesian-optimization/BayesianOptimization/blob/master/examples/constraints.ipynb + for further information. Supplying contraints is optional. + + Attributes + ---------- + optimizer : bayes_opt.BayesianOptimization + The optimizer instance of this class. + """ + + verbose: InitVar[int] = 1 + random_state: InitVar[int] = 1 + allow_duplicate_points: InitVar[bool] = True + bayes_opt_kwds: InitVar[Optional[Mapping[str, Any]]] = None + + def __post_init__( + self, verbose, random_state, allow_duplicate_points, bayes_opt_kwds + ): + """Create optimizer""" + if bayes_opt_kwds is None: + bayes_opt_kwds = {} + + if self.input.bounds is None: + raise ValueError("Input.bounds is required for this optimizer") + + self.optimizer = BayesianOptimization( + f=self._opt_func, + pbounds=self.input.bounds, + constraint=self.input.constraints, + verbose=verbose, + random_state=random_state, + allow_duplicate_points=allow_duplicate_points, + **bayes_opt_kwds, + ) + + def _target_func(self, true: pd.DataFrame, predicted: pd.DataFrame) -> Number: + """Invert the cost function because BayesianOptimization maximizes the target""" + return -self.input.cost_func(true, predicted) + + def run(self, **opt_kwargs) -> BayesianOptimizerOutput: + """Execute the optimization + + ``BayesianOptimization`` *maximizes* a target function. Therefore, this class + inverts the cost function and used that as target function. The cost function is + still minimized. + + Parameters + ---------- + init_points : int, optional + Number of initial samples taken from the parameter space. Defaults to 10^N, + where N is the number of parameters. + n_iter : int, optional + Number of iteration steps after initial sampling. Defaults to 10^N, where N + is the number of parameters. + opt_kwargs + Further keyword arguments passed to ``BayesianOptimization.maximize``. + + Returns + ------- + output : BayesianOptimizerOutput + Optimization output. :py:attr:`BayesianOptimizerOutput.p_space` stores data + on the sampled parameter space. + """ + # Retrieve parameters + num_params = len(self.input.bounds) + init_points = opt_kwargs.pop("init_points", 10**num_params) + n_iter = opt_kwargs.pop("n_iter", 10**num_params) + + # Run optimizer + self.optimizer.maximize(init_points=init_points, n_iter=n_iter, **opt_kwargs) + + # Return output + opt = self.optimizer.max + return BayesianOptimizerOutput( + params=opt["params"], + target=opt["target"], + p_space=self.optimizer.space, + ) + + +@dataclass +class BayesianOptimizerOutputEvaluator(OutputEvaluator): + """Evaluate the output of :py:class:`BayesianOptimizer`. + + Parameters + ---------- + input : Input + The input object for the optimization task. + output : BayesianOptimizerOutput + The output object returned by the Bayesian optimization task. + + Raises + ------ + TypeError + If :py:attr:`output` is not of type :py:class:`BayesianOptimizerOutput` + """ + + output: BayesianOptimizerOutput + + def __post_init__(self): + """Check output type and call base class post_init""" + if not isinstance(self.output, BayesianOptimizerOutput): + raise TypeError("'output' must be type BayesianOptimizerOutput") + + super().__post_init__() + + def plot_impf_variability( + self, + cost_func_diff: float = 0.1, + p_space_df: Optional[pd.DataFrame] = None, + plot_haz: bool = True, + plot_impf_kws: Optional[dict] = None, + plot_hist_kws: Optional[dict] = None, + ): + """Plot impact function variability with parameter combinations of + almost equal cost function values + + Args: + cost_func_diff (float, optional): Max deviation from optimal cost + function value (as fraction). Defaults to 0.1 (i.e. 10%). + p_space_df (pd.DataFrame, optional): parameter space. Defaults to None. + plot_haz (bool, optional): Whether or not to plot hazard intensity + distibution. Defaults to False. + plot_impf_kws (dict, optional): Keyword arguments for impact + function plot. Defaults to None. + plot_hist_kws (dict, optional): Keyword arguments for hazard + intensity distribution plot. Defaults to None. + """ + + # Initialize plot keyword arguments + if plot_impf_kws is None: + plot_impf_kws = {} + if plot_hist_kws is None: + plot_hist_kws = {} + + # Retrieve hazard type and parameter space + haz_type = self.input.hazard.haz_type + if p_space_df is None: + p_space_df = self.output.p_space_to_dataframe() + + # Retrieve parameters of impact functions with cost function values + # within 'cost_func_diff' % of the best estimate + params_within_range = p_space_df["Parameters"] + plot_space_label = "Parameter space" + if cost_func_diff is not None: + max_cost_func_val = p_space_df["Calibration", "Cost Function"].min() * ( + 1 + cost_func_diff + ) + params_within_range = params_within_range.loc[ + p_space_df["Calibration", "Cost Function"] <= max_cost_func_val + ] + plot_space_label = ( + f"within {int(cost_func_diff*100)} percent " f"of best fit" + ) + + # Set plot defaults + color = plot_impf_kws.pop("color", "tab:blue") + lw = plot_impf_kws.pop("lw", 2) + zorder = plot_impf_kws.pop("zorder", 3) + label = plot_impf_kws.pop("label", "best fit") + + # get number of impact functions and create a plot for each + n_impf = len(self.impf_set.get_func(haz_type=haz_type)) + axes = [] + + for impf_idx in range(n_impf): + _, ax = plt.subplots() + + # Plot best-fit impact function + best_impf = self.impf_set.get_func(haz_type=haz_type)[impf_idx] + ax.plot( + best_impf.intensity, + best_impf.mdd * best_impf.paa * 100, + color=color, + lw=lw, + zorder=zorder, + label=label, + **plot_impf_kws, + ) + + # Plot all impact functions within 'cost_func_diff' % of best estimate + for row in range(params_within_range.shape[0]): + label_temp = plot_space_label if row == 0 else None + + sel_params = params_within_range.iloc[row, :].to_dict() + temp_impf_set = self.input.impact_func_creator(**sel_params) + temp_impf = temp_impf_set.get_func(haz_type=haz_type)[impf_idx] + + ax.plot( + temp_impf.intensity, + temp_impf.mdd * temp_impf.paa * 100, + color="grey", + alpha=0.4, + label=label_temp, + ) + + # Plot hazard intensity value distributions + if plot_haz: + haz_vals = self.input.hazard.intensity[ + :, self.input.exposure.gdf[f"centr_{haz_type}"] + ] + + # Plot defaults + color_hist = plot_hist_kws.pop("color", "tab:orange") + alpha_hist = plot_hist_kws.pop("alpha", 0.3) + + ax2 = ax.twinx() + ax2.hist( + haz_vals.data, + bins=40, + color=color_hist, + alpha=alpha_hist, + label="Hazard intensity\noccurence", + ) + ax2.set(ylabel="Hazard intensity occurence (#Exposure points)") + ax.axvline( + x=haz_vals.max(), label="Maximum hazard value", color="tab:orange" + ) + ax2.legend(loc="lower right") + + ax.set( + xlabel=f"Intensity ({self.input.hazard.units})", + ylabel="Mean Damage Ratio (MDR) in %", + xlim=(min(best_impf.intensity), max(best_impf.intensity)), + ) + ax.legend() + axes.append(ax) + + if n_impf > 1: + return axes + + return ax diff --git a/doc/tutorial/climada_util_calibrate.ipynb b/doc/tutorial/climada_util_calibrate.ipynb index f6fbb42922..a0220d1203 100644 --- a/doc/tutorial/climada_util_calibrate.ipynb +++ b/doc/tutorial/climada_util_calibrate.ipynb @@ -2629,7 +2629,10 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Here we show how the variability in parameter combinations with similar cost function values (as seen in the plot of the parameter space) translate to varying impact functions. In addition, the hazard value distribution is shown. Together this provides an intuitive overview regarding the robustness of the optimization, given the chosen cost function. It does NOT provide a view of the sampling uncertainty (as e.g. bootstrapping or cross-validation) NOR of the suitability of the cost function which is chosen by the user." + "Here we show how the variability in parameter combinations with similar cost function values (as seen in the plot of the parameter space) translate to varying impact functions. In addition, the hazard value distribution is shown. Together this provides an intuitive overview regarding the robustness of the optimization, given the chosen cost function. It does NOT provide a view of the sampling uncertainty (as e.g. bootstrapping or cross-validation) NOR of the suitability of the cost function which is chosen by the user.\n", + "\n", + "This functionality is only available from the `BayesianOptimizerOutputEvaluator` tailored to Bayesian optimizer outputs.\n", + "It includes all function from `OutputEvaluator`." ] }, { @@ -2669,6 +2672,10 @@ } ], "source": [ + "from climada.util.calibrate import BayesianOptimizerOutputEvaluator\n", + "\n", + "output_eval = BayesianOptimizerOutputEvaluator(input, bayes_output)\n", + "\n", "# Plot the impact function variability\n", "output_eval.plot_impf_variability(cost_func_diff=0.1)\n", "output_eval.plot_impf_variability(\n",