Add utility module for calibrating impact functions

[peanutfun]

This issue collects ideas for an impact function calibration module

Input information

For calibrating any impact function, we need to specify four main things:

• The impact data. This is the "real world" data we want to use to calibrate the impact function. The computed ("modeled") impact should be as close to it as possible.
• The impact function. We have to specify what its overall shape is and how it is parameterized.
• The cost function. It relates the modeled impact with the impact data and returns a scalar that indicates how close they are. The task for calibrating an impact function is reduced to minimizing the cost function.
• The impact estimation task. As the impact function is varied by the optimization algorithm, we "only" need to supply the usual Hazard and Exposure.

Calibration with Bayesian Optimization

This idea is based on the Petals calib_opt module by @timschmi95, see https://github.com/CLIMADA-project/climada_petals/blob/feature/calib_opt/climada_petals/engine/calib_opt.py

The module wants us to instantiate a BayesianOptimization instance. Its main parameter is a function that has the search space parameters as input and returns the target function value. In this case, the target function value is maximized, which simply means that we have to invert the cost function value.

With regard to the impact data I would take a simplistic approach. To accomodate multiple events and multiple regions, the input data should be a pandas.DataFrame with the event_id and index and arbitrary (!) columns that can (!) represent the target impact values. The cost function then takes this data frame and the impact object of each iteration and the user can define how the cost should then be computed from both objects. In this case, we assume that the columns in the impact data are impacts per region. We compute the difference for each region and event, square all results, take the log for the sum of all regions and then sum the result of all events (as an overly complicated example for what we can do with this):

def my_cost_function(data: pd.DataFrame, impact: climada.engine.Impact) -> float:
    return np.sum(np.log10(((data - impact.impact_at_reg())**2).sum(axis=1)))

The other important input it the impact function. Its parameters are what the optimizer calibrates in order to minimize the cost function. The user has to supply a "generator" function that takes the set of calibrated parameters and returns a new impact function set. In this example, we take the TC impact function by Emanuel 2011 and intend to vary two of the three parameters:

def my_impf_set_gen(v_thresh: float, v_half: float) -> climada.entity.ImpactFuncSet:
    return climada.engine.ImpactFuncSet([
        climada.entity.ImpfTropCyclone.from_emanuel_usa(
            v_thresh=v_thresh, v_half=v_half
        )
    ])

Now we only need a utility function to plug everything together:

from bayes_opt import BayesianOptimization

def calibrate(
    hazard: Hazard,
    exposure: Exposures,
    impact_data: pd.DataFrame,
    impact_func_gen: Callable[..., ImpactFuncSet],
    cost_func: Callable[[pd.DataFrame, Impact], float],
    param_bounds: Mapping[str, Tuple[float, float]],
    optimizer_kwds: Mapping[str, Any],
    impact_calc_kwds: Mapping[str, Any],
):
    """Calibrate parameters for an impact function"""
    # Match event IDs in hazard and impact_data
    hazard = hazard.select(event_id=impact_data.index)

    # Assign the centroids
    exposure.assign_centroids(hazard)

    # Create the calibration function
    def func(**kwargs):
        impf_set = impact_func_gen(**kwargs)
        impact = ImpactCalc(exposure, impf_set, hazard).impact(
            assign_centroids=False, **impact_calc_kwds
        )
        return 1 / cost_func(impact_data, impact)

    # Create the optimizer (with some default settings)
    optimizer_opts = dict(random_state=1, allow_duplicate_points=True)
    optimizer_opts.update(**optimizer_kwds)
    optimizer = BayesianOptimization(f=func, pbounds=param_bounds, **optimizer_opts)

    return optimizer

The user can then call the optimizer like advertised by the BayesianOptimization docs. We assume we have a hazard, an exposure, and the impact data defined somwhere. Then we only need to define the bounds for the parameters we are looking for. The BayesianOptimization package automatically deduces the parameter names and initial values from these automatically.

optimizer = calibrate(
    hazard=hazard,
    exposure=exposure,
    impact_data=impact_data,
    impact_func_gen=my_impf_set_gen,
    cost_func=my_cost_function,
    # NOTE: Bad example. These parameters must be constrained because v_thresh < v_half is required!
    param_bounds={"v_half": (0, 100), "v_thresh": (0, 50)},  
    optimizer_kwds=dict(verbose=2),
    impact_calc_kwds=dict(save_mat=True),  # To be able to call impact.imact_at_reg()
)
optimizer.maximize(init_points=100, n_iter=200)
result = optimizer.max["params"]

# Use generator to receive the function we want
result_impf_set = my_impf_set_gen(result)

[timschmi95]

Sounds great to me! I think the proposed structure looks good. A few points that spring to mind:

• We need good default arguments regarding the BaysianOpitmizer options, in particular more samples the more Impf-parameters one tries to fit.
• I would suggest to also return some plots of the sampling space by default, so the user can check if it was reasonably well sampled. Not sure where to put these plotting fucntions though.
• If one wants to fit an impf-parameter with a logarithmnic dependence and test multiple orders of magnitude (e.g. the k parameter in the sigmoid function) the sampling space would have to be transformed for a well distributed sampling. This is currently not implemented, but doesn't have to be for the first version.
• In addition the display of impact functions with similar skill, a bootstrapping (of Events) could be used for more robust calibration.
• I'm not sure if it wouldn't be better to let the user input a cost_function measure (e.g. RMSE, RMSF) rather than a callable? I think if the correct matching of impact_data and imp.at_reg_id() can be guaranteed inside the function, it could be nicer than the callable. It could potentially also accept both, in case the user has a very specific cost function.

[peanutfun]

From a personal discussion yesterday:

We want to define two data classes for setting the common inputs and output of any optimization algorithm. While the input should always have exactly the same structure, the output can have additional fields that provide more information, depending on the optimization algorithm used.

The input takes the following data:

• Hazard
• Exposure
• Impact data ("true" data to compare against)
• Cost function
• Impact function generator (creates impact function from optimized parameters)
• Parameter bounds
• Parameter constraints

The output contains at least the following:

• Optimized parameters
• Cost function result at optimum
• An easy way to calculate the impact with the optimal parameter set
• (Possibly:) Additional data like the explored parameter space, depending on the optimization algorithm

The input data is to be put into a class or function that actually performs the optimization. Depending on the optimization algorithm, it transforms or wraps the data from the input to match the algorithm input parameters. We are still unclear about how this could be structured. One idea is to create an abstract base class defining the interface. This would allow easy addition of new optimization algorithms outside the Climada code base, but might lead to a confusing documentation of the classes involved. We'll have to come up with an example to see how bad or good it looks.

[thomini]

I am currently trying to retrieve impact functions based on claims we have collected. Even though I am currently not using CLIMADA for that, and it is not precisely calibration, I am sharing my thoughts here - maybe this is helpful.

Basically, I have available spatially claims data in EUR, exposure data in EUR and wind gust data in m/s. By combining claims and exposure, I get a damage ratio for each claim. That claim is associated with a storm event and a location, for which the wind gust data is available. With all that, naively I can fit an impact function to all these points (damage ratio vs. wind gust).

But here the basic question for me is about the shape and parameters of the impact function? (like mentioned above by @peanutfun)

Another usual problem when performing this kind of tasks is of course the lack of data for higher intensities (here higher wind gusts). Usually no or few data is available. Here I think calibration comes in handy, where a basic impact function as a starting point is provided.

Hope that helps!

[peanutfun]

@timschmi95 very much, thanks for the input!

I agree that we are stressing the term "calibration" a bit. It is more like an optimization or fitting procedure. But it seems to be a coined term already, so I think it's fine to stick to it.

The choice of the impact function type/shape is a critical one and (from my perspective) we don't have to tackle it in this module. Users have to decide for themselves which impact function makes sense for their use case, and we can supply some sensible suggestions, like a step function or a sigmoid function. But the choice of impact function should always be an "active" one, meaning the user has to decide explicitly.

This also relates to your last point referring to fewer data points for higher intensities. In a usual fitting procedure the uncertainty for these intensities then increases. However, the choice of impact functions greatly influences that. For example, step and sigmoid functions are monotonically increasing by definition, meaning that at least the impact will never decrease with increasing intensity -- NB: For a certain types of hazard and exposures, this might actually be what you want. Again, the choice of impact function shape is critical

[thomini]

@peanutfun thank you!

I totally agree with you, that the shape should be actively chosen by the user - finally the shape may heavily influence the high intensities and thus the results for extremes... in particular, when the impact data does not include such extremes..

Is there a tutorial / jupyter notebook where an example calibration is performed in CLIMADA?