Add option to store and load calibration results

Author: peanutfun

climada/util/calibrate/base.py

@@ -22,11 +22,13 @@
 from dataclasses import dataclass, field, InitVar
 from typing import Callable, Mapping, Optional, Tuple, Union, Any, Dict
 from numbers import Number
+from pathlib import Path
 
 import pandas as pd
 import numpy as np
 from scipy.optimize import Bounds, LinearConstraint, NonlinearConstraint
 import seaborn as sns
+import h5py
 
 from climada.hazard import Hazard
 from climada.entity import Exposures, ImpactFuncSet
@@ -185,6 +187,37 @@ class Output:
     params: Mapping[str, Number]
     target: Number
 
+    def to_hdf5(self, filepath: Union[Path, str], mode:str = "x"):
+        """Write the output into an H5 file
+
+        This stores the data as attributes because we only store single numbers, not
+        arrays
+
+        Parameters
+        ----------
+        filepath : Path or str
+            The filepath to store the data.
+        mode : str (optional)
+            The mode for opening the file. Defaults to ``x`` (Create file, fail if
+            exists).
+        """
+        with h5py.File(filepath, mode=mode) as file:
+            # Store target
+            grp = file.create_group("base")
+            grp.attrs["target"] = self.target
+
+            # Store params
+            grp_params = grp.create_group("params")
+            for p_name, p_val in self.params.items():
+                grp_params.attrs[p_name] = p_val
+
+    @classmethod
+    def from_hdf5(cls, filepath: Union[Path, str]):
+        """Create an output object from an H5 file"""
+        with h5py.File(filepath) as file:
+            target = file["base"].attrs["target"]
+            params = dict(file["base"]["params"].attrs.items())
+            return cls(params=params, target=target)
 
 @dataclass
 class OutputEvaluator:

climada/util/calibrate/bayesian_optimizer.py

@@ -24,6 +24,7 @@
 from itertools import combinations, repeat
 from collections import deque, namedtuple
 import logging
+from pathlib import Path
 
 import pandas as pd
 import numpy as np
@@ -38,6 +39,26 @@
 LOGGER = logging.getLogger(__name__)
 
 
+@dataclass
+class _FakeConstraint:
+    """Fake the behavior of the constrait for cycling the BayesianOutputOptimizer"""
+
+    results: np.ndarray
+
+    @property
+    def lb(self):
+        """Return the lower bound"""
+        return np.array([0])
+
+    def allowed(self, values):
+        """Return if the values are allowed. This only mocks the true behavior"""
+        if self.results.shape != values.shape:
+            raise ValueError("Inserting wrong constraint values")
+        return self.results
+
+
+# TODO: Add read/write method
+# TODO: Export this class
 @dataclass
 class BayesianOptimizerOutput(Output):
     """Output of a calibration with :py:class:`BayesianOptimizer`
@@ -83,6 +104,59 @@ def p_space_to_dataframe(self):
         data.index.rename("Iteration", inplace=True)
         return data
 
+    def to_hdf5(self, filepath: Union[Path, str], mode: str = "x"):
+        """Write this output to an H5 file"""
+        # Write base class information
+        super().to_hdf5(filepath=filepath, mode=mode)
+
+        # Write parameter space
+        p_space_df = self.p_space_to_dataframe()
+        p_space_df.to_hdf(filepath, mode="a", key="p_space")
+
+    @classmethod
+    def from_hdf5(cls, filepath: Union[Path, str]):
+        """Read BayesianOptimizerOutput from an H5 file
+
+        Warning
+        -------
+        This results in an object with broken :py:attr:`p_space` object. Do not further
+        modify this parameter space. This function is only intended to load the
+        parameter space again for analysis/plotting.
+        """
+        output = Output.from_hdf5(filepath)
+        p_space_df = pd.read_hdf(filepath, mode="r", key="p_space")
+        p_space_df["Calibration", "Target"] = -p_space_df[
+            "Calibration", "Cost Function"
+        ]
+
+        # Reorganize data
+        bounds = {param: (np.nan, np.nan) for param in p_space_df["Parameters"].columns}
+        constraint = None
+        if "Constraints Function" in p_space_df["Calibration"].columns:
+            constraint = _FakeConstraint(
+                p_space_df["Calibration", "Allowed"].to_numpy()
+            )
+
+        p_space = TargetSpace(
+            target_func=lambda x: x,
+            pbounds=bounds,
+            constraint=constraint,
+            allow_duplicate_points=True,
+        )
+        for _, row in p_space_df.iterrows():
+            constraint_value = (
+                None
+                if constraint is None
+                else row["Calibration", "Constraints Function"]
+            )
+            p_space.register(
+                params=row["Parameters"].to_numpy(),
+                target=row["Calibration", "Target"],
+                constraint_value=constraint_value,
+            )
+
+        return cls(params=output.params, target=output.target, p_space=p_space)
+
     def plot_p_space(
         self,
         p_space_df: Optional[pd.DataFrame] = None,
@@ -685,14 +759,18 @@ def plot_impf_variability(
                 # Plot defaults
                 color_hist = plot_hist_kws.pop("color", "tab:orange")
                 alpha_hist = plot_hist_kws.pop("alpha", 0.3)
+                bins = plot_hist_kws.pop("bins", 40)
+                label = plot_hist_kws.pop("label", "Hazard intensity\noccurence")
 
+                # Histogram plot
                 ax2 = ax.twinx()
                 ax2.hist(
                     haz_vals.data,
-                    bins=40,
+                    bins=bins,
                     color=color_hist,
                     alpha=alpha_hist,
-                    label="Hazard intensity\noccurence",
+                    label=label,
+                    **plot_hist_kws,
                 )
                 ax2.set(ylabel="Hazard intensity occurence (#Exposure points)")
                 ax.axvline(

climada/util/calibrate/test/test_base.py

@@ -20,6 +20,8 @@
 
 import unittest
 from unittest.mock import patch, create_autospec, MagicMock
+from tempfile import TemporaryDirectory
+from pathlib import Path
 
 import numpy as np
 import numpy.testing as npt
@@ -206,6 +208,20 @@ def setUp(self):
         self.optimizer = ConcreteOptimizer(self.input)
 
 
+class TestOuput(unittest.TestCase):
+    """Test the optimizer output"""
+
+    def test_cycle(self):
+        """Test if cycling an output object works"""
+        output = Output(params={"p1": 1.0, "p_2": 10}, target=2.0)
+        with TemporaryDirectory() as tmpdir:
+            outfile = Path(tmpdir, "out.h5")
+            output.to_hdf5(outfile)
+            self.assertTrue(outfile.is_file())
+            output_2 = Output.from_hdf5(outfile)
+        self.assertEqual(output.target, output_2.target)
+        self.assertDictEqual(output.params, output_2.params)
+
 class TestOutputEvaluator(unittest.TestCase):
     """Test the output evaluator"""
 

climada/util/calibrate/test/test_bayesian_optimizer.py

@@ -20,13 +20,14 @@
 
 import unittest
 from unittest.mock import patch, MagicMock
+from tempfile import TemporaryDirectory
+from pathlib import Path
 
 import numpy as np
 import numpy.testing as npt
 import pandas as pd
 from bayes_opt import BayesianOptimization, Events
 from scipy.optimize import NonlinearConstraint
-import matplotlib.pyplot as plt
 from matplotlib.axes import Axes
 
 from climada.util.calibrate import Input, BayesianOptimizer, BayesianOptimizerController
@@ -185,8 +186,8 @@ def test_improvements(self):
 class TestBayesianOptimizerOutput(unittest.TestCase):
     """Tests for the output class of BayesianOptimizer"""
 
-    def test_p_space_to_dataframe(self):
-        """"""
+    def setUp(self):
+        """Create a default output"""
         bayes_opt = BayesianOptimization(
             f=lambda x: -(x**2),
             pbounds={"x": (-10, 10)},
@@ -198,13 +199,16 @@ def test_p_space_to_dataframe(self):
         bayes_opt.probe({"x": 1.0}, lazy=False)
         bayes_opt.probe({"x": -0.9}, lazy=False)
 
-        output = BayesianOptimizerOutput(
+        self.output = BayesianOptimizerOutput(
             params=bayes_opt.max["params"],
             target=bayes_opt.max["target"],
             p_space=bayes_opt.space,
         )
-        self.assertDictEqual(output.params, {"x": 1.0})
-        self.assertEqual(output.target, -1.0)
+
+    def test_p_space_to_dataframe(self):
+        """"""
+        self.assertDictEqual(self.output.params, {"x": 1.0})
+        self.assertEqual(self.output.target, -1.0)
 
         idx = pd.MultiIndex.from_tuples(
             [
@@ -220,7 +224,19 @@ def test_p_space_to_dataframe(self):
         df["Calibration", "Constraints Function"] = df["Parameters", "x"]
         df["Calibration", "Allowed"] = [True, True, False]
         df.index.rename("Iteration", inplace=True)
-        pd.testing.assert_frame_equal(output.p_space_to_dataframe(), df)
+        pd.testing.assert_frame_equal(self.output.p_space_to_dataframe(), df)
+
+    def test_cycle(self):
+        """Check if the output can be cycled to produce the same p_space_df"""
+        with TemporaryDirectory() as tmpdir:
+            outpath = Path(tmpdir, "file.h5")
+            self.output.to_hdf5(outpath)
+            self.assertTrue(outpath.is_file())
+
+            output = BayesianOptimizerOutput.from_hdf5(outpath)
+        pd.testing.assert_frame_equal(
+            self.output.p_space_to_dataframe(), output.p_space_to_dataframe()
+        )
 
     def test_plot_p_space(self):
         """Test plotting of different parameter combinations"""